{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SciPy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The SciPy library is one of the core packages that make up the SciPy stack. It provides many user-friendly and efficient numerical routines such as routines for numerical integration and optimization.\n",
    "\n",
    "Library documentation: <a>http://www.scipy.org/scipylib/index.html</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# needed to display the graphs\n",
    "%matplotlib inline\n",
    "from pylab import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from numpy import *\n",
    "from scipy.integrate import quad, dblquad, tplquad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 0.0)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# integration\n",
    "val, abserr = quad(lambda x: exp(-x ** 2),  Inf, Inf)\n",
    "val, abserr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint, ode"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEACAYAAABfxaZOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXl8U1X6/z8nS9ukS5rue4GyyCKLrIJIAVnc+CGIIFtR\nx8EZF2ZcvoMiWJFRZ9QZBxQHVxQGGVBUFhdcKKAg4ChFsOy0BVq6L+maND2/P24TkrZpcrfkpL3v\n1yuQ5J57zunJvc957nOe8zyEUgoFBQUFha6BytcdUFBQUFDwHorQV1BQUOhCKEJfQUFBoQuhCH0F\nBQWFLoQi9BUUFBS6EIrQV1BQUOhCiBb6hJB3CSFFhJBfXRxPJ4RUEUJ+aXk9LbZNBQUFBQVhaCSo\n4z0AawB80EGZvZTSaRK0paCgoKAgAtGaPqV0P4AKN8WI2HYUFBQUFMTjDZs+BTCaEJJNCPmcENLP\nC20qKCgoKLSDFOYdd/wMIJlSWkcIuRnApwB6e6FdBQUFBYVWyC70KaUmh/dfEELWEkIiKKXljuUI\nIUoQIAUFBQUBUEo9NqHLbt4hhMQSQkjL+xEASGuBb4NSqrwoxTPPPOPzPrDyUsZCGQtlLDp+8UW0\npk8I+RDAOABRhJCLAJ4BoG0R4usA3AngD4SQJgB1AOaIbVNBQUFBQRiihT6l9G43x18H8LrYdhQU\nFBQUxKPsyGWQ9PR0X3eBGZSxuIoyFldRxkI4RIhNSA4IIZSVvigoKCj4C4QQUB4Lud5w2VRQUGiH\nFv8GBQWPkUIxVoS+goIPUZ5uFTxFKiVBsekrKCgodCEUoa+goKDQhVCEvoKCgkIXQhH6CgoKsrF+\n/XqMHTvW191wS2ZmJhYsWODrbngFRegrKCh0ebqSJ5Ui9BUUFBS6EIrQV1BQcEKlUuH8+fP2z4sW\nLcLy5csBAFlZWUhKSsI//vEPxMbGIiEhAevXr7eXLSsrw7Rp02AwGDBy5EicO3fOqe4lS5YgJSUF\nBoMBw4YNw/fff28/lpmZiVmzZmHBggUICwvDwIEDcebMGbzwwguIjY1Famoqvv76a3v59PR0PPnk\nkxg5ciQMBgOmT5+Oioqr+Zx+/PFHjB49GkajEYMHD8bevXvtxy5cuIBx48YhLCwMkydPRmlpqWTj\nxzqK0FdQYJCmpibk5uaKfjU1NYnuCyHEyfxRVFSE6upqFBQU4J133sGDDz6IqqoqAMCDDz4IvV6P\nK1eu4N1338V7773ndO6IESOQnZ2NiooKzJ07F7NmzYLZbLYf37lzJxYuXIiKigoMGTIEkyZNAgAU\nFBRg+fLlWLx4sVPfNmzYgPfeew+FhYXQaDR45JFHAACXL1/GbbfdhhUrVqCiogIvv/wyZs6cibKy\nMgDA3LlzMXz4cJSVlWH58uV4//33u4yJRxH6CgoKbnHcRKbVarFixQqo1WrcfPPNCAkJwalTp2C1\nWrFt2zasXLkSOp0O/fv3R0ZGhtO58+bNg9FohEqlwqOPPorGxkacOnXKfvzGG2/EpEmToFarceed\nd6KsrAxLly6FWq3G7NmzkZubi+rqagDcZLRw4UL069cPer0ezz33HLZs2YLm5mZs3LgRt9xyC6ZO\nnQoAuOmmmzBs2DDs2rUL+fn5+Omnn/Dcc89Bq9Vi7NixuP3227vMRjllR66CAoNoNBp069bN191o\nl8jISKhUV/VFvV6PmpoalJSUoKmpCcnJyfZjKSkpTue+/PLLePfdd1FQUABCCKqrq51MKzExMfb3\nOp0OUVFRdg1cp9MBAGpqahAWFgYAbdqyWCwoLS1FXl4etm7dih07dtiPNzU1YcKECSgoKIDRaLTX\nBwCpqam4ePGiqHHxFxShr6Cg4IRer0ddXZ39c2FhoZNwdUV0dDQ0Gg3y8/PRp08fAEB+fr79+P79\n+/HSSy/hu+++Q//+/QEAERERojRsx/rz8/Oh1WoRHR2NlJQULFiwAG+++Wabc/Ly8lBRUYG6ujro\n9Xr7d2q1WnA//AnFvKOgoODE4MGD8Z///AdWqxVffvkl9u3b59F5arUaM2bMQGZmJurr6/Hbb785\n2cpNJhM0Gg2ioqJgNpuxcuVKu6lGCJRSbNy4ETk5Oairq8OKFSswa9YsEEIwf/587NixA7t374bV\nakVDQwOysrJw+fJlpKamYtiwYXjmmWdgsVjw/fffY+fOnYL74W8oQl9BQcGJf/3rX9ixYweMRiM2\nbdqEO+64w+l4Rwuer732GmpqahAXF4d7770X9957r/3Y1KlTMXXqVPTu3RvdunWDTqdzMv+0XjBu\nry3Hz4QQLFiwAIsWLUJ8fDzMZjNWr14NAEhKSsJnn32G559/HjExMUhJScErr7yC5uZmAMCmTZtw\n6NAhREREYOXKlcjIyOA5Sv6LEk9fQcFHtMRB93U3/Jbx48djwYIFThNLZ8bV9cI3nr6i6SsoKPgt\nyqTJH0XoKygo+C1dxbdeShTzjoKCj1DMOwp8UMw7CgoKCgq8UYS+goKCQhdCEfoKCgoKXQhF6Cso\nKCh0IRShr6CgoNCFUIS+goKCrMidMtFfUh065iXwJYrQV1BQ8Gv8xVe/vTATvkAR+goKCsxii5XT\nWWBhX4Yi9BUUFNogZ8rEkydPYtKkSYiMjMQ111yDrVu3OrXzhz/8AbfccgtCQkKQlZXVpm/uUh3O\nmjUL8fHxCA8Px7hx4/Dbb7851f/HP/4Rt9xyC0JDQzF27FhcuXIFS5YsgdFoRN++fXH06FF7+W7d\nuuHFF19E//79ERERgXvvvReNjY324zt37sTgwYNhNBoxZswY/Prrr/Zjv/zyC6677jqEhYVhzpw5\naGho8HD05UUR+goKCm6RKmVibW0tJk2ahPnz56OkpASbN2/GH//4R+Tk5Njr/vDDD7F8+XLU1NRg\nzJgxbfriLtXhrbfeirNnz6KkpATXXXcd5s2b53T+1q1b8de//hWlpaUICAjAqFGjMHz4cJSXl+PO\nO+/Eo48+6lR+06ZN2L17N86dO4fTp09j1apVADihft999+Gtt95CeXk5Fi9ejGnTpsFiscBsNmP6\n9OnIyMhARUUFZs2ahY8//pgJ8w4opUy8uK4oKHQd3F7zgDQvARBC6Llz5+yfFy1aRJ9++mlKKaV7\n9uyhOp2OWq1W+/GYmBh66NAh2tTURLVaLT116pT92FNPPUVvuOEGSimlmzdvpmPHjnVq6/e//z19\n9tlnKaWUZmRk0IyMDJf9ysvLoxqNhtbV1dm/mzt3Lp0/f3675SsqKighhFZXV9v/jt///vf242vW\nrKH9+vWzfz527BgNDw+3f+7WrRtdt26d/fPnn39O09LSKKWUPvDAA3T58uVO7fXp04fu3buX7t27\nlyYkJDgdGz16dJvyfHB1vbR877GsVTJnKSiwCgP2X1cITZmYl5eHQ4cOwWg02r9ramrCwoULAXBP\nFElJSS7bdZfq0Gq1YtmyZfjoo49QUlJi72NpaSlCQ0MBOKdkDAoKapOisaamxqnN1n9LQUGB/W/5\n4IMPsGbNGvtxi8WCwsJCUEqRmJjoVE9qaqpi01dQUGCT9lImemKacEyZaMPxfUpKCsaNG4eKigr7\ny2Qy4fXXX/eoX/Hx8fZUhzby8vLsfdu0aRO2b9+Ob7/9FlVVVbhw4QIAcQuorf8WmzBPSUnBsmXL\nnP6WmpoazJ49G/Hx8bh8+bJTPY799CWK0FdQUGiDXCkTb731Vpw+fRobN26ExWKBxWLBkSNHcPLk\nSQDuhbO7VIc1NTUIDAxEREQEamtr8dRTTzmdz1f4U0qxdu1aXL58GeXl5fjrX/+K2bNnAwDuv/9+\n/Pvf/8bhw4dBKUVtbS127dqFmpoajB49GhqNBqtXr4bFYsG2bdtw5MgRXm3LhSL0FRQU2iBXysTQ\n0FDs3r0bmzdvRmJiIuLj4/Hkk0/CbDbb63WnDXeU6nDhwoVITU1FYmIiBgwYgOuvv75NisWOPrf+\n2wghmDt3LiZPnoy0tDT06tULTz/9NABg6NCheOutt/DQQw8hIiICvXr1wgcffAAA0Gq12LZtG9av\nX4/IyEhs2bIFM2fO7PDv8hai4+kTQt4FcCuAYkrptS7KrAZwM4A6AIsopb+0U4ayYO9SUPAWSjx9\n9unevTveeecdTJgwwdddYSqe/nsApro6SAi5BUBPSmkvAL8H8IYEbSooKCgoCEC00KeU7gdQ0UGR\naQDebyl7CEA4ISRWbLsKCgoKCvzxhstmIoCLDp8vAUgCUOSFthUUFBQEY/P+6Ux4y0+/tb2pXUPm\nsuXLoFVrAQDp6elIT0+XuVtXOXToEC5fvow77rjDZ25VZrMZn3zyCfr3748BAwb4pA8AcODAAZSW\nlmLatGk+60NDQwM+/fRTDBkyBH369PFZP3744QdUVlbi1ltv9VkfFBQcycrKQlZWFmpqauy7oPng\nDe+dywCSHT4ntXzXhuLkCmRmZiIzM9OrAh8ATCYTxo4di+zsbK+268jBgwdxxx132Dd/+AKb69l1\n113ntDXe2/z444+YOXOmfdONL6CUor6+Hr1790Zubq7P+qGg4Eh6ejoyMzMxdepUvPXWW7zP94bQ\n3w5gIQAQQkYBqKSUtmva2ZS3GdlXvC90CwoKkJiYiOjoaJSXl3u9fRtWqxUBAQHQ6XR2FzZvc+HC\nBfTs2RNJSUkoLCz0SR8Abiy0Wi3UarXPIi0ePXoUQ4YMQa9evdoEDVNQ8DWOO6J5nSe2YULIhwAO\nAOhDCLlICLmXELKYELIYACilnwM4Twg5C2AdgD+6qitiz2O4f8f9sDZbxXaLFydPnsQ111zj1TY7\nYuDAgU7R+rzJhQsX0L17d5+0bYNSajexDRw4EMeOHfNJPyorKxEZGemTthUUOqK2thbBwcGCzpXC\ne+duSmkCpTSAUppMKX2XUrqOUrrOocxDlNKelNJBlNKfXdVVefBhBJJgrDm8xlUR2bAJGbVajaam\nJq+3f+XKFcTHxwMADAaDIFudFLRez/CFH/mlS5fs8VoiIyNRUdGRc5h8OP7tQrUqBQU5yM7OxuDB\ngwWdy9SVPFL9E+arXsGqfauQW5nrtXYdBV2PHj18smJ/5swZ9OrVq90+eRPHdpOTk3Hp0iWv9+Hs\n2bNIS0vzersdERoaiurqal93g0n+85//YMqUKS6P79+/X5In6W7duuHbb7/1qKzcKRpdkZ+fj9DQ\nUNmVJbPZjMDAQEHnMiX0B0acR+FeIx4f/Tge2PmAT7TMpKQknwg6q9UKjca3QU9b2867d+/uM5c1\nXwemMpvNCAgIsH/u16+fTxe2WWbevHn46quv7J9bJ2AZO3asPbaOGFhJN9gRKSkpMJlMTPeTKaHf\nv1sJjvxE8dj1j+FKzRVs+nWT7G2WlJQgKirK/pmlH8vbk15eXh66detm/6zRaGC1end9hRXOnz+P\nHj162D+3jjqpwOHq+lDCS7ALU0J/xLhg/HQ+AhqVFm/d/hYe2/0YSmpLZG3z9OnTTmYVX9H6JomP\nj8eVK1e82of8/HykpqZ6tU1PCAsL87pppbCw0L7G0hW5ePEiZsyYgZiYGERFReHhhx8GwJlNxowZ\ng0cffRRRUVHIzMx0MqXceOONAIBBgwYhNDQUW7duRVZWllNMeld1nzt3DhMmTEBUVBSio6Mxf/58\nj9e23KVoXLJkCVJSUmAwGDBs2DB8//339mOZmZmYNWsWFixYgLCwMAwcOBBnzpzBCy+8gNjYWKSm\npuLrr7+2l09PT8eTTz6JkSNHwmAwYPr06fZ1p9zcXKhUKvtTc3p6OlasWIEbbrgBYWFhmDJlCsrK\nyux1ffDBB0hNTUVUVBRWrVrlsQlLzKTKlNDvNy4NtMmKy5eB4YnDMe/aefjzV3+WtU2z2YygoCBZ\n2/CE1k8YPXr08LqbIKW0zYKltzW29trr0aOHk7nAW/1g6anPm1itVtx2223o3r078vLycPnyZdx9\n993244cPH0ZaWhqKi4uxbNkyp3NtIZiPHTsGk8mEWbNmua17zpw59uPLli1DYWEhcnJycPHiRWRm\nZnrU545SNALAiBEjkJ2djYqKCsydOxezZs1ycoveuXMnFi5ciIqKCgwZMgSTJk0CwLlzL1++HIsX\nL3Zqb8OGDXjvvfdQWFgIjUaDRx55xGXfPvzwQ6xfvx7FxcUwm814+eWXAQC//fYbHnzwQXz44Yco\nLCxEVVUVCgoK3F53NTU1CAkJ8Whc2oMpoa/umYbB6mP4pSUG58rxK/HjpR/x6clPfdsxmWlqaoJa\nrXb6LiAgABaLxUc9uoq3BV9lZaVTViUAMBqNqKys9Go/WIAQaV58OXz4MAoLC/HSSy9Bp9MhMDAQ\no0ePth9PSEjAgw8+CJVKxVthaq9uWx7ctLQ0TJw4EVqtFlFRUfjzn/+MvXv3uq3TarVi27ZtWLly\nJXQ6Hfr374+MjAwnBWLevHkwGo1QqVR49NFH0djYiFOnTtmP33jjjZg0aRLUajXuvPNOlJWVYenS\npVCr1Zg9ezZyc3PtT5uEECxcuBD9+vWDXq/Hc889hy1btriMgHnPPfegZ8+eCAoKwl133WVPvP7R\nRx9h2rRpGD16NLRaLVauXOnR/Zafn++UjYwvTAl9dOuGgZafceJX7tEoOCAY709/H3/Y9QcU1xZ7\nrRtardarAregoAAJCQlea48v3tT2L1686GQK8BUsaPlSJcnly8WLF5GamurSTVXM79NR3UVFRZgz\nZw6SkpJgMBiwYMECJ1OIK9ylaASAl19+Gf369UN4eDiMRiOqqqpQWlpqP946ZWJUVJT9GrClZnRM\no9i6LYvF4lSfI3FxcU512+opKChwSg2p0+k82hdSXFyM6Ohot+VcwZbQDwrCNfoLOP5Tvf2rMSlj\nsGjQIvx+x++9Jny87cHDiqBrj5iYGBQXe2/CLS8vR0REhNfa40N4eLjP9gx4k+TkZOTn57tcpBUz\nIXZU91NPPQW1Wo3jx4+jqqoKGzZs8Gg3trsUjfv378dLL72ErVu3orKyEhUVFTAYDJKmULQ9nfAh\nISHBSc7U19d7NMm1Z4blA1tCH0Cv2BIcP+b8Y2SmZ+JC5QV8kP2BV/qQmJjoVaHf2j2QJVJTU5GX\nl+fVNlnQstsjKSmpTd7TzsjIkSMRHx+PpUuXoq6uDg0NDThw4IDH58fGxrpcjxoxYoTLumtqahAc\nHIywsDBcvnwZL730kkftuUvRaDKZoNFoEBUVBbPZjJUrV4pyDKCUYuPGjcjJyUFdXR1WrFiBWbNm\nubxuXU0uM2fOxI4dO3Dw4EGYzWZkZmZ6RbFlTugnp9bhVH4QHDfFBmoCseGODXji6yeQVymdAHI1\nwFqt1ie7cn1Je+sKABASEoLa2lqv9YMFV7/6+vp2bdWRkZEuH+E7EyqVCjt27MDZs2eRkpKC5ORk\nbNmyBYDr9IKO32VmZiIjIwNGoxEfffSR03G1Wu2y7meeeQY///wzDAYDbr/9dsycOdNjBaCjFI1T\np07F1KlT0bt3b3Tr1g06nc7J/OMuZWLrz4QQLFiwAIsWLUJ8fDzMZjNWr17t8bm2z/3798eaNWsw\nZ84cJCQkIDQ0FDExMYI3XXkMpZSJF9cVSi9kZNA0YynNyaFteHH/i3Tsu2OpxWppe1AApaWlNDs7\nu91j3333nSRteIKrtrzZh9zcXHru3Dmf98NVW3v27PFaH86cOUPz8/PbPSblWNiueQX/Ij09nb7z\nzjuS12symahGo6G5ubntHrddL62vwZbvPZa1zGn62h490Fd/HsePtz32xJgnEKgJxHN7n5OkLbGr\n4HITEhICk8nklbYuXbrktKjkCAvmFpVK5bWNYgUFBV3aR1/BPVSiJ9IdO3agrq4OtbW1ePzxxzFw\n4EDZ98owJ/TD+vbFNTjertBXERU23LEBb/38FvZc2CO6rcrKShgMBtH1yIU3N2hZLBaX6wpSXeCe\n4GqC8eaCMgshMRTYRipFaPv27UhMTERiYiLOnTuHzZs3d1ieSrB/hDmhH9KnD/o3/tyu0AeAuJA4\nrJ++Hgs+WSDJbl0WtFhXxMbGoqio62SVpFdNfW1ITEz0aXIZBQUbe/bscVozEMNbb72FiooKVFZW\n4uuvv3YbHaCioqLNPha+MCf0SWIiBtUeQkfh5CenTcb8gfOR8WkGmqlvEmxIhdVqbXcBFfD+fgFX\neGtirK6uRlhYWLvHQkNDvWbq6giWlQSFzk9RURFiY2NF1cGc0EdEBK6xnEBeHkVDg+tiz41/DpUN\nlfjb93+TpRveurlLS0t5+/d6G2+Zd4qLizu8oL1pZlJQYJHWASKFwJ7QJwSINqBHkgUOu6TboFVr\nsWXWFqw5vAZfnf3KdUGBeGsRVYqZW240Go1XnjiKi4uddkayCCHEZ+kbFRSkWG9iT+gDaIyKQv+k\nKpw40XG5pLAkbL5zMxZ+uhDnK6QNyBUXF+eVHLHl5eWibXRyExMTg5ISeaOdAh0vJrNCV/HVV+i8\nsCv0I6+4FfoAcGPqjVg2dhlm/HcG6izSxTuPjY31ircIFbml2hvEx8d7ZQJ0Z77xlsmto3a8NQEq\nKMgFk9LGHBmJ/voLHgl9AHh4xMO4NvZa3L/jfsnsvqwsonoDd2PWlVIFWiwWlwvrAKfpexIfRUE6\nWmfiYpHWcfRZhkmh3xgVhf7kN4+FPiEE625bh5ySHLx84GWPzmloaJB/u7Of4M4NrCt5rJSWlnYY\nwVCtVvvFja2g4AomhX5zdDR61B/DpUtAfb378gCg1+rx2ZzP8K9D/8K2nG1uy7vzFGEFbywcFhUV\nOYV/7cr4w8J6Z6GrxbeSAikUMCaFvq5bN5DCfKSlAXzyKScbkvHZnM+weOdiHLl8pMOyntzcLGi4\n4eHhHqeME0pJSYlHcbzlxt14e8OLqLKyEuHh4bK24Q/k5OQgPT0dRqMRAwYMwI4dOwAAhw4dQnx8\nvJNJ8JNPPsGgQYMAAM3NzXjxxRfRs2dPREVFYfbs2W1SCb777rtITU3FTTfd1G7bL730EhISEpCU\nlIR3333X6diuXbswZMgQGAwGpKSk4Nlnn7Ufs9W/fv16pKSkIDIyEv/+979x5MgRDBw4EEaj0Z6a\nEbia+vHhhx9GeHg4+vbti++++85+vKqqCvfdd5+9L8uXL7crYM3NzXj88ccRHR2NtLQ07Nq1S8xw\nexUmhX5Yr15ovnIF/fvDYxOPjaEJQ/H27W9j+n+ndxiRU2zKMW8RFRUlu7eIJ25gLEyA3hgLgI2/\n1ZdYLBbcfvvtmDp1KkpKSrBmzRrMmzcPZ86cwciRIxEcHOyUx3XTpk2YN28eAGDNmjXYvn079u3b\nh8LCQhiNRjz44INO9e/btw8nT57EV1+1dbX+8ssv8corr+Cbb77B6dOn8c033zgdDwkJwcaNG1FV\nVYVdu3bhjTfewGeffeZU5vDhwzh79iw2b96MJUuW4Pnnn8d3332HEydOYMuWLfaUjrayPXv2RFlZ\nGZ599lnMmDHDnqVt0aJFCAgIwLlz5/DLL79g9+7dePvttwEAb775Jnbt2oWjR4/ip59+skcTlRtJ\n1iz5RGeT8wWHiIPWggLaGB5On32W0qVL3YSmc8E/D/6TDlg7gFbWV7Z73JNoid6I7OiuH/X19fTA\ngQM+7QOl3hkLd22UlpbSo0ePytoHb44F3ETZRCYkefFl3759NC4uzum7u+++m2ZmZlJKKX366afp\nvffeSymltLq6mgYHB9ujkvbt25d+++239vMKCgqoVqulVquVXrhwgRJC6IULF1y2fc8999Ann3zS\n/vn06dOUEOIyAuySJUvon//8Z0optddfUFBgPx4ZGUm3bNli/zxz5kz66quvUkopfe+992hCQoJT\nfSNGjKAbNmygV65coYGBgbS+vt5+bNOmTXT8+PGUUkrHjx9P161bZz+2e/duSgihVqvV5d8mFgD0\n8OHD7X5PechaJqNKqaKjoTGZ0P8aK97f6NqToiOWjFyCc+XnMP2/0/HFvC8QpPF98nMhBAUFobGx\n0dfdkB1X8fwdMRqN+LWj+BydDPqMb3YgFxQUtMnklpqaak8gc/fdd2PMmDF44403sG3bNgwdOtRe\nPjc3F3fccYeTG7JGo3GKIdVRlrjCwkIMHz7c/rl1FNxDhw5h6dKlOHHiBMxmMxobG3HXXXc5lXE0\n2+p0ujafHfNDJCYmtvk7CwoKkJ+fD4vF4hRttbm52d6fwsLCDtMzyoUUZlgmzTvQaNAUGor+8eW8\nzTs2CCH4183/QlxIHGZ/NBtNzewtGikeRFcpKytze0GrVColFIMXSEhIwMWLF53GOi8vzx56u1+/\nfkhNTcUXX3yBTZs2Ye7cufZyKSkp+PLLL1FRUWF/1dXVOQnPjswg8fHxLtMeAsDcuXMxffp0XLp0\nCZWVlXjggQdEOTq0zoSWl5eHxMREJCcnIzAwEGVlZfa/o6qqyq50uOunXHReoQ/AbDSiZ3AhCgqA\nOoF7rlREhfenvw+z1Yzfbf+doOBscgqZsrIyj+JoyC3oWLBhs5wbt6sxatQo6PV6/P3vf4fFYkFW\nVhZ27tyJOXPm2MvMnTsXr776Kvbv349Zs2bZv3/ggQfw1FNP2YVgSUkJtm/f7nHbd911F9avX29P\nRei4UAtwa3FGoxEBAQE4fPgwNm3axPv6dbyfiouLsXr1algsFmzduhUnT57ELbfcgri4OEyePBmP\nPvooTCYTmpubce7cOft6wF133YXVq1fj8uXLqKiowIsvvsirD0JxFZCQD8wKfUt4ODRlRejVC8jJ\nEV5PgDoAH836CKfLTuOJ3U/wEqAGg0HW+DueaLesILfnTFlZmV8J/c7sq6/VarFjxw588cUXiI6O\nxkMPPYQNGzagd+/e9jJ333039u3bh4kTJzr9bkuWLMG0adMwefJkhIWF4frrr8fhw4ftx90J6KlT\np+JPf/oTJkyYgN69e2PixIlO56xduxYrVqxAWFgYnnvuOcyePdvpfE8mAMcyI0eOxJkzZxAdHY3l\ny5fj448/tu9Z+eCDD2A2m9GvXz9ERERg1qxZ9vwW999/P6ZMmYJBgwZh2LBhvFI7ikGSNvgsAMj5\nQqtFrSstzjdiAAAgAElEQVQTJlC6YQOdM4fS998XsuzhTFldGR2wdgB9Zs8zlFLPFuzy8/Pp2bNn\nxTfugr1799Kmpia35eROV+hJ/Tk5ObSwsNCnfeBTTs5+HD16lJaWlopuq/U1r+Bd3nvvPXrDDTf4\nuhse4+p6gb+nS7RhNhqBoiJBbpvtEaGLwDcLvsGWE1uwcu9Kj86R20WwubnZ7eIlK3jLXdIdLJii\noqKiUF5e7utuKCgIglmhb4mIAIqLJRP6ABAbEos9GXuw+fhmbMjb4La8TqdDQ0dB/b0EC4LOaDQy\nIeiojOsbntbNylgoiIMQwsS95W2YFfqN4eGSavo2YkNi8V3Gd/i2+FuPEqzLKWRYoKPMXY7IHXOG\nhZuvrq4OwcHBbsvpdDrUCfUuUGCGjIwMp41aXQVmhb6lxbyTlgYUFQE1NdLVbVAb8PqI17Hp+Cas\n2LOCecEuZ/9YCTvAwm/gqQdRV9UQFToHzAp9xMaiuagIajXQu7c4D57WlJWVoXdCb+xdtBc7Tu/A\nki+X+CTXrqeCTs6ga/7kQQTI+0Tgb2OhoCAE0UKfEDKVEHKSEHKGEPKXdo6nE0KqCCG/tLye9qRe\nfbduoC3uUVKbeGw3d0xwDLIysnD0ylFkfJoBi5XN+PkGg0G2ePZ8BB0L2q1arZYtOmNVVRUMBoNH\nZVl4MlFQEIIooU8IUQN4DcBUAP0A3E0I6dtO0b2U0iEtr1We1B3WsydIaSlAqeRC39GkYQgy4Mv5\nX6K8vhwztsxAvcXDWM5exGg02iMVSk1DQwOCgjwLUcGCoDMajfaAWFLT3Nzs9SxmNlOR8lJe7l5S\nIfYKHwHgLKU0l1JqAbAZwP9rpxzvHhvj49Gs1QLV1Rg0CMjOFtlTB1q7Suq1enw6+1OEBYZh4gcT\nUVJ7NR2elIMtlIiIiE7vLeLpOHcmzxlXftR79uzxyN/6zJkz9nAJUr887QPfsnxf3333nUflsrOz\nUVpa6tM+yDkWFotFskVnsUI/EcBFh8+XWr5zhAIYTQjJJoR8Tgjp50nFWq0WlvBwoLgY110H/Pwz\nQGVUNLVqLTbcsQETuk/AqHdG4WQpF8ifytmoh4SGhsq2M5iFSY0PrEyAco1bc3MzExMgC9d9XV0d\ndDqdR2XlfBpmAXfZ7fggVuh7cmX8DCCZUjoIwBoAn3pauTk8HCgpQXw8oNUCFy+6P0cMKqLCqgmr\n8PTYpzFu/Thk5WbJ1panrpIAJ2BYuAnlorGxEQEBAR6V1el0qPc0nZofwmddITw8XDZTFwuUl5d7\nvN7EijIgF1I6GYgNrXwZgGOc1GRw2r4dSqnJ4f0XhJC1hJAISmmbXygzM9P+Pj09HQNaNH0Adm3f\nGxFM7xlyD1LDUzH7o9m4J+kepCNd8jYqKys9vrnlhIXJhG+wNRaeTuQaNz43t5x7J1gY47Kysg7D\nMDui1+tl2zvBwliUl5ejR48eAICsrCxkZWUJrkus0P8JQC9CSDcABQBmA7jbsQAhJBZAMaWUEkJG\nACDtCXzAWegDQGErof/LL8D06SJ77CETuk9AVkYWJr83GTWf1+AfU/6BALVn2qgnVFZWMiHoWLmg\nO0pG3pUoLy/HgAEDPC4v1+TDgjJQWVmJa6+91qOyLFzHcmI2m+1Pw+np6UhPT7cfax2J1B2izDuU\n0iYADwH4CsBvAP5LKc0hhCwmhCxuKXYngF8JIUcBvApgTvu1tcVm3gGuavpS4OkF0je6LzaO24jz\nZecx8YOJKDQVStMB8N8UxcJNqNVqZYm06Y9hleUyufGxY3d2fOFN1R4s3Hssee+AUvoFpbQPpbQn\npfSFlu/WUUrXtbx/nVI6gFI6mFI6mlL6o6d1mw2GNuYdb5MSk4LXbngNk3tMxvC3huPAxQOS1Gsy\nmZjI0cvngg4PD5dlsaypqcltjl5HWLgJQ0JCnDIwSQWllNcNzoKGq1KpYLVaJa+X798m13XBwhhL\nie+n0Q6wGI12TT8lBaivB1r2a3mNiIgIVFVWYfm45Xjz9jdxx3/vwCsHXhG9g9cfb+7O7iHBB7k8\nZ1j4nQF+/ZBz74Sv4TuRBAUFMe9owLTQNzvY9AkBhg4F/vc/7/YhNDTUvhv2ll634NDvDuHjnI9x\n66ZbUVRT5OZs1/ijFiOX0PdHQceKt4gc1wVfsworeyfkuI4aGhp4mdsiIiKYV4yYFvoWB5s+AIwa\nBRw86N0+tM7L2i28G/Yu2ouh8UMxZN0Q7D6327sdkhC+AqMzh5rmm6/YYDCgqqpK8n6wQHV1Na+0\nfP4g6ITCd+2NlQmwI5gW+o6aPgCMHi1e6PPZ/OIKrVqLVRNWYeOMjbj3s3vx6FePos7Cz12Mr8CV\nQ9DV19czsWjIgo2ej388wE6SdjmuC76CrjOHmhYi9Fk3dTEt9JvCw0FLS4EWX+RRo4AjRwAx8bZM\nJpMkyYUBzq3z6ANHUVhTiMH/HizZIm97yCFgWAmrzBdlLK7CwliwsnlQrrHgowxoNBrZAgJKBdNC\nPyQiAggOBlpmTqMRSE4Gfv1VeJ1S39xR+ih8OPNDvHjTi5i5ZSYe++oxWYK2yXExVVZWSra125vI\nIWSqqqp4XxcsCDo54Gve6cw0NjZ6HJBQTqS81pgW+kajEU0taRNtXH+9OBOPXBrdjL4z8OsffsVl\n02UMXjfYbQgHvo/lcmy5r6ioYGJXMN+x0Ov1kntI+Kt/vBzukqz4x/srrCsDTP+yERERaDQYnBZz\nR48GDoiwolRVVcmmxUTpo7D5zs34201/w8JPFmLBJwtEefg4IofnjNls5rV4CbBxQRsMBlnspqx4\nEfHBYDDIFoyPDyyMXVBQEBOOBqzDtNAPCwtDQ2hom8XcH34QXiefQGc2+Aq66ddMx28P/ob4kHgM\neGMA1h5ZC2uzszbGt045XARZuFGFIMdTj5DJjIXxk2sC5AsrykBn9KiSwvnEEaaFvkql4jR9B6Hf\npw+3SSs313v9EDLgIQEh+Pukv2NPxh5sPr4ZI98eif15+wX3oTN7SPAVGJ315rbFTudDZ420abFY\neO3SBrixYOG6kFoZqKmpQWhoqGT1MS30gba++oQAEycC33zjw07xYEDMAOxdtBePXv8o5n8yHzP+\nOwNnys7wvjCkzp4DsKGdWa1W3vbjoKAgNDY2ytQjz5F6/Orq6hAcHMzrnLCwMNlSafJB6mtTyMI6\nK089UiP1OqR/CH0HTR8AJk0Cvv7ae30Qe3MTQjD32rk4+eBJjEgcgevfuR5rz69FWV2ZRD30X0wm\nk6DFZKkFLgumGiE3t1qtliXuja8RMhZ6vV6WeEi+pusJfYf4OzYmTgS+/dbuvu8VpBAyOq0OS29Y\niqy7sqAOUKPPa33wbNazqGrwzSOpvwo6VpB6/Px5LKRGyFjI8TQsBDnciaV0PmFe6Jtb2fQBzlc/\nKgo4epR/fUIuCqldBLWNWrwy4RX8+Lsfcb7yPHqu6YlV+1ahutH3j+nukNpFUKigY+HmlhqhiXWk\nHgsh9Ukt6GpqanibuuToBwtYrVbe6xsdwb7Qb2XTtzFpkvfs+lIvltk2RfWM6In3p7+PH+79ASdL\nT6Ln6p54fv/zqGzwjl1SyA0SFhYmqYtgVVWVpItUQhEyFhqNRtL8AhaLhbcLLcCOoOuMJjchsN5v\n/xD6rTR9AJgyBdi1S8KGTp0CFiwAHnkEKC11OiS1V0B9fb3TLr/ekb2xccZG7F20F7+V/Ia01Wl4\nfPfjuFTtlHmSiZtb6gmwubmZtwstwMZYdGYvIr7o9XomfORZELiEENnSWEoB80LfGh4OlJcDrUwK\nN90EZGe3+xDAn8JCYMIEoH9/wGzmZhQHc47UXgGuYun3je6LjTM24pfFv6CZNmPgGwOx6NNFOF58\nXLK2HRFyg0g9AbJwkwJsjIVQWBjDzug5I1SxkPppWGqYF/pUrQZsgt+BoCBg8mRg+3YJGlm2DJg3\nD1i6FHjjDaBbN+D55x3a8q6LYIohBf+Y8g+ce+Qc+kT2waQNkzBpwyTsL92PpmbfBnMKCwtjQtCx\nQGf0kRe6EagzjgXfWPo2WJ8AmRf6AIDo6HZNPHfcAXzyici6Cwq4SpYt4z4TArz6KvD660CZb10q\njTojnhz7JHKX5OKewfdg66Wt6P6v7nhu73Oi8/UKffxUq9VMPLpKqd0K2QgEcCkTWdbohCA0Ci0r\nTz1SItTJgPUJkHmhTykFYmLatePccguwfz/gaUiadgNJbdgA3Hkn4Og1kZzMzSivvSai567hK7AC\nNYGYe+1crBmyBjvv3olL1ZfQb20/3LnlTuw6vUuQ9i/1Lj+hsGCb5xtL3wYrIYWlRKig64xxb4SO\nBetPw8wLfQCgLjR9gwGYOhX47389q6dNyFhKgfXrgXvuaVv4sceAdevswftZubkHxQ3CutvXIXdJ\nLib1mIRV+1ch+Z/JeGL3E7xs//4aVtmGlL+HGP94FuzpUqLsFbiKUBdajUbD9IY55oV+SEgImtrZ\noGUjI4OT257Q5oI+cwYwmbh4za3p14/Lxi7D1l+hAsvRR94QZMDiYYtx8L6DyMrIglatxdSNUzHs\nzWH458F/4mLVxQ7rUm7uqwjZ8t9ZEbMRiAXFSMo+iImlz7IywLzQNxgMqG8VadORyZOBvDzg5En3\ndbURdF9+yT0quPqBFi3yfEbxAq7irPSJ6oPnJz6PvD/l4YWJL+BEyQkMWTcEo98Z7XICqKysZCJR\nhtCbIzAwULLFdZPJhJCQEEHnsiDoAOn6IdSFVsEZVq6L9mBe6IeHh6NWr3cp9DUa4N57gbVr3dfV\nxrxjE/qumDWLKyNxwg6hgs7dYplapcaktEl4e9rbKHysECvGrcDx4uMYvG4wRr8zGi8feBknS0+C\nUir5Lj9v01ldR4UQHBzMRARWFsaws8YikhLmhb7BYEB1UFCHDvl//COwcaM9q6JLnBZyLRZuFXji\nRNcnREYCw4YBu3czcUHzcQXTqrWY2nMq3vl/76DwsUIsv3E5zpWfw6QNk9BrTS+8fu51fHP+G5it\nZpl7LQ9SusWJ0cr87bqQExa0W1aijkqF1LH0AT8Q+oGBgagLDnap6QNAYiJw883AW2/xqDg7G+je\nnUu82xF33AFs28ajYveI2fQh5IIOUAfg5l43443b3kD+n/Lx8V0fw6A1YPme5Yh5KQYz/jsDa4+s\nxanSU16/cYW2x7pbHF+ExNK30dncJYW60AKd77qQw8uOeaEPtB9pszX/93/AP/4BeBxZ9ccfgVGj\n3JebPh3YuRNEoqTkYoSqFI+uhBAMihuE+SnzcfC+gzjz8BnM6DsDRwqOYNKGSUj+ZzIyPs3AB9kf\n4HL1ZVFtucPVzmRP0Ov1kpk0xGhSUk2S9fX10Ov1gs5lRdOXSiMVs7DOygQo1VjI4XDhF0bd9iJt\ntmbQIGDcOGD1auDJJz2o9McfudAL7khKAnr2RFh2Nhf7QST19fWCogfKRXRwNOYPnI/5A+eDUoqz\n5Wfx7YVvsf3Udvz5qz8jQheBMcljuFfKGFwTdY1kbYvRYlgwq0iJmJs7MDAQZrN/munaQ8xYhIaG\ndirzTkVFBZKTkyWt0y+EviU0FKiu5nzmO3jsW7kSGDOGW9iNjW173EkrO3gQeOopzzpw++2IOnCA\n890XiVDfXxtSCbv26iGEoFdkL/SK7IUHhj2AZtqME8Un8MPFH7Avfx9e+P4FVDVWobeuN6app2FU\n0ihcF38dDEHC/p7O4DYqpUYXGRkpSV3+TmVlJa65RphywcqOcamoqqrCgAEDJK3TL4Q+1GogIoKL\nfhkX57JY797cPqvHHuMWdl1SXMyFWPD0wpo2DVFr16LJYoFGq+XX91b4082tIipcG3stro29Fg8M\newAAUGgqxJtfvoni2mI8vedpZF/JRkJoAoYmDMXQ+KEYljAM18Vfh7BA9+6glZWVSEhIENw/qUwr\nLCxAVlVVoUePHr7uBhOmLqGx9KWGhbGQw4XWP4Q+cDX+TgdCHwCeeQYYMAD4/HMuTEO7HD4MjBgB\neJqb9dproSYENYcPI3zMGH79bgUrN7dQ4kPjMSV5CgYPHoygoCBYm604WXoSPxX8hP8V/g/bcrbh\nWNExJIQmYFDcIAyIHoABMdwrLSINGtXVS06MRscKtjC6fPP8tkbMRqDOSGcz37GEXwh9QojL+Dut\nCQ4G3n8fuOsu4MgRLoxOG7KzgSFD+HQAjZMno/nTTzn7kQg6w81tWziMi4uDWqVG/5j+6B/THxmD\nMwAATc1NOFl6Er8W/YrjxcfxwbEPcLz4OApNhegT1QcDYgbg2phrYSmyILoqGj2MPRCgDuDdDxYE\ng82jyt/NVFKg0WjQ1NTk1/s/pMIWl4mFa7Q1fvHrUEpdRtpsjxtvBP70J2DGDGDPHqDNZstjx4Bp\n03j1gUyfjoDly4GXXuJ1HquIefy0ucXFuXjq0qg0du3ekRpzDXJKcnC8+Dh+Lf4VBwoPYP2H63Gx\n6iISwxLRK6IXekX0Qu/I3tzaQkQvpIanOj0dyIGYG9PmLaII/atJZXxtvpTCtCK2DlsAOiGhmeXG\nL4Q+AI81fRt/+Qtw9iwn+D/9FHDyhsvOBp5+mlfz+ilTYJ0/H7hyxa2JSTC5uVxkz/37uUXr664D\nfvc7YORIexEW7M8GgwH5+fm8zwsJCMHwxOEYnjgcAJAVmIX09HRYrBZcqLyAM2VncKb8DHJKc7D9\n9HacLjuNKzVXEB8Sj9TwVKQaUtEtvBtSDalIDU/FlboraGxqRKCGf4pBqQgPD0d+fj5SU1N91gep\nELsRyKYM+FroS4FYgW1TBhShLwYemj7AhdP597+5hd2bbnJItlJXxwXr6dOHV/OqoCCUDhuGmF27\ngPvu43WuR7z/PrcCfd993IYDrRbYt48L+zx2LBffX6KomGIfO3U6naTb/rVqLXpH9kbvyN5tjpmt\nZlysuoi8qjzkVeYhryoP+/P3Y+OvG3Hqyinc9/N9iNRFIiksCQmhCfZXfEi80+dIfSRUpK3d3Wq1\nirLHsx5Glw9iNwKFh4fj7NmzEvZIGFKYVMR6ljmaQFnDf4R+TAzw88+8TtFoOFm6bBkwbBjFn/5k\nxPjQ3ziBH8Dfhlw6ejRiduwQJfTb1dRfew145RVOw+/b9+r3I0ZwMSaeegoYOhT47DPB7TpSU1Mj\nOMAYIJ0t3ZOnlgB1ANIi0pAWkdbmWFZWFsbeOBYFpoI2r+8vfo8CUwEKTYUoMBWgurEacSFxiA+N\nR0xwDKL10YjWRyNYFQzUAjVnahClj+K+D45GsDbYo7+zM8V6ESvoWEkqI8XTcGVlJSIiIgSfHx4e\njry8PNH9kOPJXrTQJ4RMBfAqADWAtymlf2unzGoANwOoA7CIUvoLzzZ4m3dsqFTACy8AgwfX4qGH\n+uP4NefxbM9xSORdE1A2YgQnoOvrAake23bt4lIzHjjApWlsjV7PZfIaMQK46SYYli0Dxo8X1aTQ\npCEsolapkWxIRrKh4w0sDU0NuFJzBQWmApTUlqCkrgSldaU4W3gW5Q3lOHj4oP37ktoSUFBE66MR\npY+CUWdEeFA4jEHc/47vjTojzlefR1xJnP2YTsveI70nVFZWIikpSfD5KpWKCXu6FFRWVqJ79+6C\nz2c5qYwooU8IUQN4DcBNAC4DOEII2U4pzXEocwuAnpTSXoSQkQDeAOBB/AOHTmo0sISHQ8vDvNOa\nkSNL8c03BJvvy8e1n7+I2zOAxYu5SAyePt03hYVxXj/ffgvcdpugfjhpjwUFnP3ps8/aF/iOzJ0L\nREVhwF13gfboASKwfYC7oGPb273mZ/ARDkGaIHQL74Zu4d2cvj969Ci6devWRsOtNdeitK4UpXWl\nqGyoREVDBfd/Pfd/oanQ/l1eUR5ev/S6/TgAhAaGIjQgFCEBIQgNbPk/oNX/DmXyivPQcLYBoQGh\n0Gv10Gl10Gl00Gl13GeNDmqVvCGPq6qq0L9/f1F1sCCwpXgSNZvNorzsWPTasSFW0x8B4CylNBcA\nCCGbAfw/ADkOZaYBeB8AKKWHCCHhhJBYSmmRp42Eh4fDZDIhQoCmb6OqqgopKSl4wfh3PPYuwfuX\nb8LvfselWrz5Zi6Ew9Ch3H6tDj3Opk3jFghECF0AXNauxYuBBx5oP4lLe0yejNOvvIL+997LxZuY\nM0dQ01VVVejdu639HJcucUlj9u3jVsGLirgF5bAwLqpd797A8OHcU4cEN7fYG0MKF8E24bZbCA4I\nRnBAMFLD3S/QZmVxC9I26i31qDHXwGQ2wdRosr+vMdfA1Ghyel9oKoTJbML5kvM4+ONBmBpNqLPU\nob6pHvWWeqf3GpXGaRKwTQx6rR71pnokXknkjmv0CNIEIVATiAB1AALUAQhUc+9t37X+HKAOQE5p\nDlSXVR2W0ag09ld7ayRSwILAZGHykguxQj8RgGOGjksARnpQJgmAx0LfYDCgsqEBESI0fXtGoJMn\nETUyDY9159ZNz53jNnJ99RVnZbl4kUuYlZrKveLiuPVToxG4dCkKSL0LQX99AIG/b0agToWgICAw\nkHup1dwCskp19eX4mRDAauXkJdm0CcjPBz7+mNffob3hBpRu3ozoBQuAqipu4uCJ2WxGgOOaxuHD\nwN/+BmRlcave48dzCWTi4rgZsLoauHwZOHGCizj6l79gTG0tkJ4O3HAD9xoyhFt8FkNtLff0Y3uV\nlwMNDdzLbL460EFBgF6PlLIy1BUWIiwxEQgNdX4FeubR0+HGKquVa7uxkXuZzVd/ULWae6lUUFdX\nc2PU8r1OrYYuKALR+ijXCXpa0XriaA2lFGar2WkScJwYfsr+Ccndk9FEmlBnqUNDUwPMVjMamxph\ntppRa6lFRUPF1e+arx5rtHL/F5UWYWv5VqfvWpexNlthabagqbkJKqJymgQ0Kg1oE4XuF53Td1qV\ntk05Vy+tWovS4lK8V/keVEQFFVRQq9Tce6KCmnDv2/vO9n1ebh72791vn5jclW/v+5PFJ1H6WykI\nCAghICBQEZX9veP/KqJq8x0BwdHyo0Au2j3maV3na88jqiiqw/J8ISJjic8EMJVSen/L5/kARlJK\nH3YoswPAi5TSH1o+fwPg/yilP7eqi7rqi8lkwqmcHAwbM4YTDAIWYffs2YPxw4dzawMmE3fDtkNt\nLefck5vL/V9UxMXpr6gAzp4thUYThcbD2WhM6YlGdTAaG6/Khebmqy9KXb2naG6+Kgj4KzUOY0Sp\nQwV86nSowxanhJCWE12f7FQvbb7aB/vvRrjT7XW5ap46/U9AnZ8cbOd6Uoerzx123naO/R8FnlBQ\ngDQDqibXL7Wl1XdW+3vapnyrsqQZIJQ7hzQDxPZ/s8N3bb+n7X7f7LJ8x99buT6AOvzf3HKLUOdj\npLntd6CgHRyz/40e1tW2PLj/V58DpdRjSSJW078MwHEFLRmcJt9RmaSW79qQmZlpf5+enm7XfEJC\nQlBTV8e5bZaUcKYGnhBCgNOngZ49XQp8gNvR268f92pNVtZxrk9LP+S02uee492PrKy9SN+/H/TX\n46CbPczo7kBDQyOOHj2KUaNGAZcLODPTNddwSdzDwtxbXSwW/JaZiX67dnGC8IknuAxhbrT01vVm\nZe1z1kwrK4FDh4AffuCC2V26BBQWArU13HGi4uInRUVxv2P37jivUqH7+AlAWhqX2yA6mtcsWF1d\njbNnz+K6665re7CxkZvcTSagpoZ7mUzcHxIYyCkOAQE4kp2N4aNHw+mRzfbSaDzqz969ezFu3Lj2\nB81qdf1qbrbX/8OBAxgzZkz7E17r79r5v6CgAPX19Ujr2bPjsh0c2//99xg7dqxnbQIA0QBwVsDc\nPbG4o6mpCT/++CNuuOEGz05o54I/ffo0goODkShATthw+Zt60L6Nffv24cYbbxTcB5PJhNOnT2Po\n0KHOfdu3F3v37bV1AKuwile9YoX+TwB6EUK6ASgAMBvA3a3KbAfwEIDNhJBRACpd2fMdhb4jti3N\ndl99oT/myZOeB1nriNtv51wpBQj9gLIy4NVXQX76CUJMonp9EBob67nF5+RE4McD3PbjQddyJppZ\ns9qf1EpLgbffBt54A0mRkVC/9CIwZYqQRw0AgFpNndc+osKBW6dwr/ZweirhnnguZ2WhtwhPJKMx\nBA0NpvbnK20gEBIIxEd1WIe5sQbagX07LOMOjYa6mDMJuFus49usoaEBgclx0CbGCO5DZJAaOTk5\n0EYKz3us0mmgDRG30U0bqII2ULit31Rbg4gog6g6ouMiUVpaCm2Q8IVvbZBa1PkAoA5UQasTLmLr\nymoRFRfZpo6bpkzETVOuZvxb9bwXhT6ltIkQ8hCAr8C5bL5DKc0hhCxuOb6OUvo5IeQWQshZALUA\n7hHcoEC3TTsihb7d/DRqFKfF5ua697pxoKmpCd3Xr+diP4twB3NCp+O0/L17gaVLuWwyt90G9O/P\nae+XLnH+/z//DMycCWzbhmyTSZQ2BghYbGtVXoq8AiqVSnQYXRYW7KQI4xAcHIxajzMIsYsUY2Ew\nGERvEmPhupDLy060nz6l9AsAX7T6bl2rzw+JbYcQwntXbqs+AKdO8Y650y5qNXDrrcCOHcDDD7sv\n30LNwYNcXP4NG8T3oTXjxnFmlWPHgO++4xZdGxuBhARuxXr8eM52BXALtj5GbF4BqZDCU0SsgJBi\nLOxPwyLrEIsUYyE2Cm1gYCAaGxtF1SEFYsfCpZedSPxnRy4gjab/f/8n+HSnyHnTpnEbtTwV+pQi\n4MknYXrkEUTIGZxr4EDuxThSxWgRK6hY0OiqqqokyY7EwliI7YPJZJI8J6y/0sbLTiL8IkcuwD/S\nZhuam4EzZ3jH3HHE6RH65ps5rdrTx8gvvoDq0iVoHhL90MMEKpVKVPgBVjR9KRCrZdfV1TEZmMsX\nSJGbAGDD15+FPrSH3wh9MaEYACCouJjzHhERc8Yp6XJQELeb9o033J9osQCPPYYzixcjRKKgab4m\nLEZCJ0QAABjhSURBVCxMVJwVsTsepUKKGzMkJES0PZ0FASFFH2xJZRTEI9c14TdCHwAn9AVq+sEX\nL4r23AkPD0dFRcXVL/7wBy6im7sbfu1aIDERZaNGdRotxhZFUOFqSGF/RwrzTmhoKGpqaiTojf+j\nVqvR1NTk6260wb+Evgjzji4/X7TQbyPouncHJkzgQiK44sIFzrXztdcEu0eyiNNTjwBYsKUD0vSD\nFaEv5m+R6vdgZSzEYLFYJMn+JfYekQv/EvoizDv6/HxR9nwACAgIgMVicf5y1SouLHJ7k1FTExeG\n+YknpNkfwBCdIY68VOnsOsNTj1RZnjqD0K+urpZkvUnsWMilGPmX0Beo6VNKoZfAvNMuvXsD99/P\n2fcdbZmUAo8/zrl3PvaY9O0KRCotRq1Wi7LdsmCiEps0xIZWq22rDPgZYmPp2+gMyoBUY8HqBOhf\nQj8sjAt6VV/P67T6+noEX7okidBvV1itXMlt8b/nHi5IT1kZl+Zw/37gv/91E7aTP2I0ACWf61Wk\nurkBNiYxMUg1FhqNxu+Tykg1FmKTynT5hVyVSgVrc/PV+Ds8qL54EZraWuHhG9yh1XKhOrVaICmJ\nC8+pUnGboERk33GFmItBSldJFmzIYuhMbqOAuOuiMyXWEYvYzHI2pNgwJwd+I/QNBgP32CjArl9/\n9CisaWmeZ0sRQkgIF9umqooLPvbWW1yIXxkQ4xUgpXbLAmJCMSiC7ioNDQ1MuNCygFRrPaziN0Lf\nbh8T4LbZdOIEVO2FzZQDjUZyc05rxNgKpdJixCLVTRUaGir4EbqpqQlasTkAWmBBo2PVRdBXsPCb\niKHLL+TaBZ2AxVzt+fNQi0wDxxJihL6UWgwL2hCri2W+gFUXQV+g0+mYzVHra/xG6IeGhqK6ulqQ\neUcKd02WMBqNiqBrQewE2Jnw9wlQyt/D38eisbERgR5mf+OL3wh9lUolOP5OsAQbs1giODiYiV2P\nLCzksuIiKPSpx2q1Qt1BUh8++Lugq62tlcz06O9jIed6k98IfTuxscCVK56Xb2pCUEEB0KuXJM0L\nFVZSajGsegV4ipRajJj9AiyYp1wlZheCWBdBXyOloPP3DXNyOlz4n9BPTuYSg3jKhQswR0YCer18\nffKAuro60UlDWEOo0FT2ClxF2StwFSnHIigoiImY+kKR8x7xT6Gfn+95+ZMnUZeSImkXhGjZUs/c\nLNzgQp82WHEblfJpyd/HggUUF9qr1NbWQi+ToupXQp8QclXT9/SRXmKhLzSKYGfUboX6yHe2TVE2\nhAh+qZOGsGD2E9oHuZKG+CtdfkeuHb2e2wjlqQePxEJfqOcMK4JOygvJYDBwHlU8kVOL4YOUY6HX\n6wW5CEqVNIQllJj6VxF6jck5efvn1ZaSAly86FlZiYW+UK8AqbUYFjQ6oWNBCGHCPCUl/u4tIiVK\nTP2raDQa5oLx+afQ99SuTymQkyOp0Pd3rwApaZNUxkewMIEoQv8qrIyFEMVI6icUoRvm5Lym/VPo\ne6rpt5iAzBKaVViJIsiCoAsLCxNk3pH6KUVIfZRSSfvBijIg5LqQKty2DTFPgL5GqnDbNliZAB3x\nT6HvqaZ/8iS3KYuBi4kFpBa29g1zfkh9fb2k6wpBQUF+u+1faq8ZoRMgC9eS1A4XitCXCk81/ZMn\n0dy7t6RaDCsIuUHkCLYmpB8saHSseFN1xrFg5WlYCFI7XAjdPa8s5LZgjyLIQ9OvT01lwmuGBVgR\ndCwgh3+8EAHOgnbLyl4BViZAKeUFC39Ta/xK6NtdBFNSPBP6v/0GU1ISExe01AQGBvI2J8jhNsqC\noBMSloIVF1oWMJlMTITbZmECtFgskoXbFoOykNuC3T6WkMClJXSXNjE7G6UJCZLf3CyYNIxGI2+v\ngM664zE4OBi1tbW8zmElaYjU14XQmPosaqRSwMJEwhr+KfTVaqBHD+DsWdeFi4uBhgaUBQdLuhov\nFKkvPiELRBaLRfIdjyxMgELGQg4hJ9SLSEqEbJhjRTBK/ZsIianfWSc/R/xK6DtFEezdGzh92nXh\nY8eAgQPRTGmn2/EICPORZ+XmlhoWPSR8BStjwYLwFDIWLNwjVqtVMe/YcHIRdCf0s7OBQYNk6QcL\nF3RQUBDq3Zm3vAALNn17/mQf9gHgPxZy5GIV4i7JwvVsNpslt6WzMgHyRW6HC78S+k54IvQHDvRe\nfzpALgHD92Zl4eZuamqSLGmIDa1W65e5YWtrayUPt23PMMcDFrTbiooKGI1GSev01/SRcoyFI51X\n6B87Jpumzxc5bm4hsHBzd9bFZCFUVFQgIiJC0jqFbJhjQRkoLy+XfCz8dcOcHGPhiH8L/VOnuPg6\nrWls5CaE/v1luaD5RhGU4+b2VyoqKhAZGSl5vSxMaHwpLy+XVaPzJ+TWbv2JmpoaWZVE/xX6MTGc\nF8/ly22P/fwzF35Br5dFGPCNIijXBc2CRsc3imBnFnSs7BVgQXPnixJL3xllIbc9CAGGDgX+97+2\nxw4cAEaPlq1pvgtErGwEkmMC5DsWrMTSl4PAwEBeKfpoJ/Uss+GPT19dAf++4oYOBX76qe33jAn9\n5uZmyRcvATY0Or5jwUosfTn64K8Lh3KgxNS/Cmsx9QULfUJIBCHka0LIaULIbkJIuz5GhJBcQsgx\nQsgvhJDDwrvK4aQ9tKfpUyq70A8LC/PLm1sOQSc0k5jUsDCRhIeHo7y83Nfd4KVhy6WN+6O7ZGNj\nIwIDAyWvl68yIPe1LEbTXwrga0ppbwDftnxuDwognVI6hFI6QkR7bbFp+o4Xbk4OEBAApKZK2pQj\nGo2GCRdBvje3HDe4Xq/nHQLB11gsFlmevCIiIpgQ+nyQOkevDaPRyMRY8Lnm5Vp7Y20CFCP0pwF4\nv+X9+wCmd1BWnqkrKQkIDQV+/fXqd59/DtxyC0AIrFZrp7aZ8qGurk4WjwAhwc58TWVlpSw3d1BQ\nEC+bPgvI5VkmZMOcHPDRmuUaC1YmQBtiJGIspbSo5X0RgFgX5SiAbwghPxFC7hfRHoBWPyIhwM03\nc4LexrZtwO23A1BCCTvS2V3i+Ew8cvpBs2Bm4tMHubyp1Gq13yVHl2ssWHsa7jC7CCHkawBx7Rxa\n5viBUkoJIa7uujGU0kJCSDSArwkhJyml+9srmJmZaX+fnp6O9PT0NmVsIYXtERJnzAAeeQT4y184\n005uLjB5MgB5/eNZubk93cpfUVGB+Ph4L/SKfSoqKpAiYd5kf6aqqgphYWGy1M3CEyCfPsjlNiq1\nrMjKykJWVpbg8zsU+pTSSa6OEUKKCCFxlNIrhJB4AMUu6ihs+b+EEPIJgBEA3Ap9V0RGRqKsrAyJ\niYncF+PHcxr/J58A777LTQAtmbLKy8vRt29ft3X6K7ZHaE+eZiorKzv1WOh0OtTX10On07kt29DQ\n4FE5f8W2edBT06ZiAvUvWivEzz77LK/zxfza2wFktLzPAPBp6wKEED0hJLTlfTCAyQB+bV2ODzah\n79AI8MYbwD33ALW1wKOP2g/JvbPN17QZiw5oamqSLW0kC089ERERHo9FZ4ePPV1ObZyF6yIgIABm\ns9nX3WBiLGyIEfovAphECDkNYELLZxBCEgghu1rKxAHYTwg5CuAQgJ2U0t1iOtzuSvgNN3BJVfbs\n4Tx3HGBpsKWGj9DvzOMAcGPh6WIZC2YHOeEzFnLi6TjL+XvwuUdYoLGxUfadyYJVP0ppOYCb2vm+\nAMCtLe/PAxgsuHft4HKBiOFHVLkuar1ej7q6Op/2gU/dcoQStmEwGHD8+HFZ6uaDp2NhsVhke/KK\niIjAiRMnkJaW5rYsC8qAyWSSbV0hMjISRUVFfrOe5Q2HC3YlZSeBj22VLyzcsHyora2VLRerv3mL\nyOlkwIoy4On1KWc8JqPR6HGyITnvJ0/H2RvBGTu10GfhR+wKoYT53NwsRBtlYbLs7G6jgOf3iJxj\nodFoYLVaZalbDrwRkLBTC305tRhPM1d1BUHnaajpzhxhky+sjAULE4Sc5h1/Q4481q3p1EJfTjxd\nICotLUVUVJQXeuQ7PPUWkdMnnA9yKgOeeovIFeeFL3Kbdzypv7m5mYnJh4UFfm/0oVMLfTkvJE+F\nfn19facNJWwjMjISpaWlbst19lDCgH/G35GLsLAw3qkbOyueKgPemPw69x0oI54GUZJ75mZBQ2LF\nLY4FTY2VsfAEucfLn/ZOyLmPBfBcGVA0fYZhxVvEk4tE7gtJr9d7lIuUhQnKarXKEmHTBmvBtTpC\nriB8NljZL+DJdVdWViZLGk8bLCkDnVboNzU1yXpz+xNyukraYGHy8QS5F9a1Wq1HYbdZmABLSkoQ\nHR0tW/0hISEeJVKReyw8ue7kXntjyeznl0LfE28RuWduf0Lum9ufUMbiKnKPhacLuSwoA2VlZbIq\nA6zk4AD8VOh7Yk8vLS1l4uaWW4vxxHW0K3gQAZ4JGUUZuEptbW2njk3FB6vVKqtN31OUhVwXeGIf\nY8UPWm6ioqLcjoU3kpGzYK7wRBnwxs3NwlhotVpmvEVY6IM7ZYCFcQCUhVyXeCL0Wbm55f4RPV0g\nYuWilpOYmBgUF7cb4Zs55L4uYmJiUFJSImsbUiBXGk9HDAaDW9dRFkxMcnsQ2fBLoa/T6TzyFpEb\nFi4UVtLSuRsLb4xVdHS0W0HHwuQnZ7A1GzExMSgqKnJf0MfU1NTIkqPXkdjYWCbGwt215y3To18K\nfX/BGx5EKpWKCddRd1RXV8u+G9dTzxm5cTfBeWONhZWNUe4EnTcW1j1RBljAW2tvfiv0WdDYgI5v\ncG9EzGMFFm5ugI2nL3d4YyxYuT/cUVJSIrug02q1sFgssrYhBd6K0+W3Qp+FmzssLAwmk8nl8aKi\nIsTGusoX37nQaDQd3lhdyVXSnRdRV1IG3CH3BjFPYWGSVGz6fkBcXByuXLni8nhpaWmXcQ90t4jq\nlMy+k+POi6grxCCyERgYyMT6mzu8oUS621/krYmna1x5MuFO0HWlmzsuLs4vFsu8cXPHxsa6vS68\nAQtPw/7gUWW1Wr1yn7KyK9dvJRILF7S7hUMWHhnlTFHoSGhoKBMLhx1dF956fI6KimJ+4dBkMske\nmgNwPwGyQHFxsVfMsKx4Efmt0O+I+vr6LmNKcIfc28ttuJtYWJgAvbXGwooXUUcUFhYiISFB9nY8\njb/jSwoLC72SQ5eVoGt+K/SDgoJc2gq99SOygkqlcpkS7sqVK4iLi/Nyj9rCwpNZQUFBl7ouOqKo\nqAgxMTG+7obX6Oj681ZyH1Yi8/qt0O/oUYkVQectOtIgutJiMsBt3HOVFNwb0UZZoiNloCvFmvEE\nVvrhDfxa6LvynDGbzUykovMWHXkRsbKY7K2bKjk5GZcuXfJKW0Lx1liwYkN2hTef/lhfc/JmKHjf\nSwOBBAcHu9ToWDAleBPWE3d4yzsCcD8B+hpvrjclJyfj4sWLXmlLCBUVFV4LipiSksLEWLia8L25\np8dvhT7r1NXVQafTeaUtlUrFhEBzhbe8IwB27Kau8NYCKsC+YuTNsegoFIM3TTuuxt2b65CK0BeJ\nq3j2XW0xGXBtQy4oKPDazc0KrqJHFhYWdqn1JsC1oPPmelNHu6RZmAC9EXjOhiL0RZKSkoL8/Pw2\n33dFQRcfH4/CwsI235tM/7+9s4uNMivj+O+fLgXaBRehWSlT0mlpC4ilEJBGY0zMXhCjq3ux6kbj\nxhiv/FiNMYoXxku9MGpivFB3N2vULaaKWbOKru6CJiZGwkeB0tIvylCyQP2iNHT46OPFvDMOw7SF\naTnnHd7zS0jmfZl5z9Mz5zzznOfjnClnAzouzFWIc/PmTWpra53JEYcA5VzB9bjEm1wyX3DdmQxe\nW18kcRjQc/mQXRUCxYm4+5BdEueA8uzsrNO5E+e+yGazTn+EGxsbyxpGLqlqpV9uWTY7O+vUeojz\n1sbZbNZpFtOKFSvIZrPO2rsfXOVi54lzcN11SvNchlEc3CqZTIaNGzc6ay+VSnk3jKpa6ZfzpyfR\nlw65AGbpsjGTydDU1ORJov8Th8k9NjZGS0uLs/bmsqTj0Bfj4+M0Nzc7a6/c2PRBue/EtRt25cqV\nZQ0jl+OiqpV+uWXj+fPnnf5yz4Xryd3Y2MjFixfvuJfEuALktrwuPU3MZXpgnjgo+GXLlt2lZHzs\neBqHvii3y+Xs7Kyz/Pi4UNVKv1zgMKkDupw/PQ4D2sUZqKW0tbUxNDTktM1ylFqWrn3pAK2trYyM\njDhtsxylf7eLIyNLSaVSdxmJPuZuaZvXrl1zeqZAVSv9uORklw5oH9kqy5cvj4U/vXRAuzoNqJhV\nq1bNe7iNLyYmJpyvvOJalTs8PExra6vTNpubmzl37pzTNu+FgYEBNm/e7Ky9qlb6EA8re+3atUxO\nThauBwcH6ejo8CiRP+rq6pieni5cJ7kvampq7thtc3R01Lmii0OGW57iueoj9haX2EIpLs6PLqbq\nlX4caG9v5+zZs4Xrqakpp19intIJ7uMHccuWLZw5c6Zwnc1mnVUmF1P8t/syDNLpNKOjo4Xr27dv\ne3e3+ZIhlUoxMTFRuJbk/QdpZmbGyx5dpUaiaypW+pKelnRa0m1JO+d5315JA5KGJH210vbmeX5h\nUk9PT3tRMKXbPPtSMsXtTk5OejmHdfXq1XdsbBWHlZjrtLw85XzIPij+DnysNgBaWlruiC3EYY74\nWoW2t7czODjovN08i7H0TwJPAX+Z6w2SaoAfAHuBrcAzkrYsos276OjoKFiWx48fp6uraykff88s\n5SA+dOhQRZ9rampifHwcgFOnTrFt27Ylk+l+WEoru9K+aGhoKJzYNDIyQjqdXpQclbDUK69K+6I4\nmymTyXjJ6Cp2rSxFcL/SvihO8/a17XhxPYuP4H7FSt/MBszs7AJveycwbGbnzOwm0AN8qNI2y7F+\n/fpC4Uc2m/V2YlZ+xXH9+vVFrzYqHdDFmRpxyNy5cuUK69atW9QzKu2LrVu30t/fD7gv2Csmr9wy\nmQypVGpRz6q0Lzo7O+nr6ytc+3arDA0N0dbWtqhnVNoX27dv58SJE4D/fgA/xtmDngkbgOI8wgvR\nvSUln8HjM5Mnv+I4cuQIO3bs8CJDPg/ZzLwGrPJL+b6+Pjo7O73IkA+i+s7uSqfTDA8PMzg4SHt7\nuxcZamtrmZmZ4datW14V3Zo1a5icnPRaS1NfX8/U1BTZbNbrNin19fVcvXqVy5cv09DQ4LTteZW+\npNcknSzz74P3+HwnjrtNmzbR09PjTcEAbNiwgbGxMW7cuOH1fN7m5mZ6enrYtWuXNxnS6TT9/f3U\n1NR4VTKpVIr9+/ezZ88ebzK0tLRw+vRpamtrvfZFY2Mjvb29dHd3e5Ohq6uLw4cPU1dX500GyB1c\nf+DAAa99sXv3bg4ePOjlJDct1rcm6Q3gy2Z2tMz/dQPfNLO90fU+YNbMvl3mvf4jfoFAIFCFmNk9\nWxRLtb6Zq8EjQJukZuAi8FHgmXJvvB+hA4FAIFAZi0nZfEpSBugGXpX0++h+o6RXAczsFvA54A9A\nP7DfzM7M9cxAIBAIPFgW7d4JBAKBQPXgvSL3QRdvVQuSmiS9ERW8nZL0Bd8y+UZSjaRjkn7rWxaf\nSHpMUq+kM5L6o1hZIpG0L5ojJyX9QpL7klpPSHpB0iVJJ4vuvTVKuDkr6Y+SHlvoOV6VvovirSri\nJvAlM3s7OZfZZxPcF3meI+cWTPpy9PvA78xsC9AJJNJFGsUGPwPsNLN3ADXAx3zK5JgXyenKYr4G\nvGZm7cCfo+t58W3pP/DirWrBzN40s+PR62vkJnbyNsOPkJQC3g/8hLkTBR56JL0FeI+ZvQC5OJmZ\n/XeBjz2sXCVnHNVJegSoAybm/8jDg5n9Ffh3ye0ngZei1y8BH17oOb6VvpPirWojsmh2AH/3K4lX\nvgt8BfC/d7Zf0sAVSS9KOirpx5L8Jrp7wsz+BXwHOE8uG/A/ZvYnv1J553Ezy++dfQl4fKEP+Fb6\nSV+234WkR4Fe4LnI4k8ckj4AXDazYyTYyo94BNgJ/NDMdgLT3MMS/mFEUivwRaCZ3Cr4UUkf9ypU\njLBcVs6COtW30p8Aig9xbSJn7ScSScuAXwE/M7Pf+JbHI+8CnpQ0BrwMvE/STz3L5IsLwAUz+0d0\n3UvuRyCJ7AL+Zmb/jNLBf01urCSZS5LeBiBpPXB5oQ/4VvqF4i1JteSKt17xLJMXlKvRfx7oN7Pv\n+ZbHJ2b2dTNrMrM0uUDd62b2Sd9y+cDM3gQykvIb9zwBnPYokk8GgG5JK6P58gS5QH+SeQV4Nnr9\nLLCgsehvxyFyQSlJ+eKtGuD5BBdvvRv4BNAn6Vh0b5+ZHfQoU1xIuhvw88DPI8NoBPiUZ3m8YGYn\nohXfEXKxnqPAj/xK5Q5JLwPvBdZFhbHfAL4F/FLSp4FzwEcWfE4ozgoEAoHk4Nu9EwgEAgGHBKUf\nCAQCCSIo/UAgEEgQQekHAoFAgghKPxAIBBJEUPqBQCCQIILSDwQCgQQRlH4gEAgkiP8BYYh8Ti/Q\nduUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xaa65fd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# differential equation\n",
    "def dy(y, t, zeta, w0):\n",
    "    x, p = y[0], y[1]\n",
    "    \n",
    "    dx = p\n",
    "    dp = -2 * zeta * w0 * p - w0**2 * x\n",
    "\n",
    "    return [dx, dp]\n",
    "\n",
    "# initial state\n",
    "y0 = [1.0, 0.0]\n",
    "\n",
    "# time coodinate to solve the ODE for\n",
    "t = linspace(0, 10, 1000)\n",
    "w0 = 2*pi*1.0\n",
    "\n",
    "# solve the ODE problem for three different values of the damping ratio\n",
    "y1 = odeint(dy, y0, t, args=(0.0, w0)) # undamped\n",
    "y2 = odeint(dy, y0, t, args=(0.2, w0)) # under damped\n",
    "y3 = odeint(dy, y0, t, args=(1.0, w0)) # critial damping\n",
    "y4 = odeint(dy, y0, t, args=(5.0, w0)) # over damped\n",
    "\n",
    "fig, ax = subplots()\n",
    "ax.plot(t, y1[:,0], 'k', label=\"undamped\", linewidth=0.25)\n",
    "ax.plot(t, y2[:,0], 'r', label=\"under damped\")\n",
    "ax.plot(t, y3[:,0], 'b', label=r\"critical damping\")\n",
    "ax.plot(t, y4[:,0], 'g', label=\"over damped\")\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy.fftpack import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhkAAADICAYAAABF5/MoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHORJREFUeJzt3X+QXGWd7/H3J5PJLwIEDCQRsiuuICC4icoPF10bhK0U\nxSKuW4oWGr2sy+IuUtRdFdjadaK1V/FeUUtXdpUfFbn3cpcSNguCSkTa1bouAW4CkZAFXAKJMpPw\nIyG/M8l87x/P6UzPpKf7dM/0dPfM51V1qk+fc7r7yUnS85nv85znKCIwMzMzG2tTWt0AMzMzm5gc\nMszMzKwpHDLMzMysKRwyzMzMrCkcMszMzKwpHDLMzMysKXKFDEldklZLujd73iNpU7ZttaQlzW2m\nmZmZdZqpOY+7GlgHHJ49D+DGiLixKa0yMzOzjlezkiHpeOBC4GZApc1l62ZmZmaHyNNd8jXgM8BA\n2bYArpL0uKRbJM1pSuvMzMysY1XtLpF0EbA5IlZLKpTtugn4Qrb+ReCrwOUVXu85y83MzCaQiMjd\nk1GrkvEHwMWSngPuAM6T9L2I2BwZUjfKmVUa42UMls9//vMtb8NEWnw+fT7befH59Lls16VeVUNG\nRFwfEQsj4gTgUuCnEfExSQvKDns/sLbuTzYzM7MJLe/VJZAGepZizFck/X72/DngirFumJmZmXW2\n3CEjIopAMVv/aJPaYyMoFAqtbsKE4vM5tnw+x5bP59jxuWwtNdLHkvvNpWjm+5uZmdn4kUSM4cBP\nM7MhPvEJ+O1vR95/113wne+MX3vMrH25kmFmdZHge9+Dj47QaXrssbBlC/i/vtnE40qGmTXdwMDI\n+6b4W8XMMv46MLO6VatSOGSYWYm/Dsysbq5kmFke/jows7o5ZJhZHv46MLO6OWSYWR7+OjCzulUL\nGV1d49cOM2tvDhlmVjcP/DSzPHJ9HUjqkrRa0r3Z86MlrZT0tKQHJM1pbjPNrJ1Uq2Qo9xX0ZjbR\n5f2d42pgHYM3SLsWWBkRJwEPZs/NbIIrVTAOHBj5GFcyzKyk5teBpOOBC4GbSXdiBbgYWJ6tLwcu\naUrrzKyt7N+fHvfuHfkYhwwzK8nzdfA14DNAeYF0XkT0Zet9wLyxbpiZtZ/+/vS4b9/IxzhkmFlJ\n1Vu9S7oI2BwRqyUVKh0TESFpxGFgPT09B9cLhYJvu2vWwUohw5UMs8mhWCxSLBYbfn3VG6RJ+m/A\nR4H9wAzgCOBu4AygEBG9khYAD0XEyRVe7xukmU0gL70ExxwDn/0s3HBD5WMWL4Y1a3yDNLOJaExv\nkBYR10fEwog4AbgU+GlEfBS4B1iaHbYUWNFog82sc7i7xMzqUe/XQel3ky8DF0h6Gjgve25mE1ye\ngZ++hNXMSqqOySgXET8DfpatvwKc36xGmVl7ciXDzOrhrwMzy60UMqrNk+FKhpmVOGSYWW4OGWZW\nD4cMM8utFDJKYzMqKV1V4qtLzMwhw8xyK4WLapWM0n1Nqh1jZpODQ4aZ5ZanklHaV+0YM5scHDLM\nLLc8YzJKx5QezWzycsgws9xcyTCzejhkmFlueSoZDhlmVuKQYWa57d8P06a5kmFm+dQMGZJmSHpY\n0hpJ6yR9KdveI2mTpNXZsqT5zTWzVurvh5kz81UyPCbDzGpOKx4ReySdGxG7JE0FfiHpXaT7mNwY\nETc2vZVm1hb6+2HGjNqVjJkzXckws5zdJRGxK1udBnQBr2bPPbef2SRSChm1Khm1goiZTQ65Qoak\nKZLWAH3AQxHxZLbrKkmPS7pF0pymtdLM2kLeSsaMGe4uMbP8lYyBiFgEHA/8oaQCcBNwArAIeBH4\narMaaWbtodQVUmueDHeXmBnUcat3gIjYJuk+4B0RUSxtl3QzcG+l1/T09BxcLxQKFAqFRtppZm2g\nVMnYvXvkYzwmw2ziKBaLFIvFhl+vqHEXI0lzgf0RsVXSTODHwDLgyYjozY65BjgjIj4y7LVR6/3N\nrHN84xuwYgW89BKsXVv5mK4uWLwYvvUtOPvs8W2fmTWXJCIi93jMPJWMBcBySVNI3Su3R8SDkr4n\naRHpKpPngCsaarGZdYxaAz8HBtIyfborGWaW7xLWtcDbKmz/WFNaZGZtq9Z4iwMHYOpU6O52yDAz\nz/hpZnUoXTkyUiVj//4UMqZO9dUlZuaQYWZ1qFXJKA8ZrmSYmUOGmeVWa0xGf7+7S8xskEOGmeVW\nazIud5eYWTmHDDPLrVYlY//+VMVwd4mZgUOGmdUh75gMd5eYGThkmFkd6rm6xCHDzBwyzCw3j8kw\ns3o4ZJhZbnnGZLi7xMxKHDLMLLdaYzJKl7C6u8TMoEbIkDRD0sOS1khaJ+lL2fajJa2U9LSkByTN\nGZ/mmlkrlSoZAwNQ6d6H7i4xs3JVQ0ZE7AHOjYhFwFuBcyW9C7gWWBkRJwEPZs/NbIIrXaI6ZUrl\nLpPSfneXmBnk6C6JiF3Z6jSgC3gVuBhYnm1fDlzSlNaZWVvp7x+cB2OkkOHuEjMrqRkyJE2RtAbo\nAx6KiCeBeRHRlx3SB8xrYhvNrE2Uxlx0dVUOEe4uMbNyeW71PgAsknQk8GNJ5w7bH5Iq9M4mPT09\nB9cLhQKFQqHhxppZa+WtZLi7xGxiKBaLFIvFhl9fM2SURMQ2SfcBbwf6JM2PiF5JC4DNI72uPGSY\nWWcrhYw8lYxduw7db2adZXhxYNmyZXW9vtbVJXNLV45ImglcAKwG7gGWZoctBVbU9alm1pHK701S\na0yGu0vMrFYlYwGwXNIUUiC5PSIelLQauFPS5cAG4IPNbaaZtYNaYzJ8q3czK1c1ZETEWuBtFba/\nApzfrEaZWXvKMybDd2E1sxLP+GlmudUzJsPdJWbmkGFmuXmeDDOrh0OGmeVW6g6pVcnwmAwzA4cM\nM6tD+Q3QXMkws1ocMswsN4/JMLN6OGSYWW61xmT4ElYzK+eQYWa51VPJcMgwM4cMM8stz4yfpf3u\nLjEzhwwzyyUi/11Y3V1iZuCQYWY5HTgAU6akxVeXmFkeNUOGpIWSHpL0pKRfSfp0tr1H0iZJq7Nl\nSfOba2atUhqPAb66xMzyyXOr937gmohYI2k28JiklUAAN0bEjU1toZm1hfKQUa2SMWuWu0vMLKkZ\nMiKiF+jN1ndIego4LtutJrbNzNpIqUoBte/C6kqGmUGdYzIkvQFYDPx7tukqSY9LukXSnDFum5m1\nkb17Yfr0tF5rTMa0abBv3/i2z8zaT57uEgCyrpLvA1dnFY2bgC9ku78IfBW4fPjrenp6Dq4XCgUK\nhcIommtmrbJvXwoPUH1MRne3Q4bZRFEsFikWiw2/PlfIkNQN3AX8z4hYARARm8v23wzcW+m15SHD\nzDpXPZWM6dPT8WbW2YYXB5YtW1bX6/NcXSLgFmBdRHy9bPuCssPeD6yt65PNrKOUh4xaV5c4ZJgZ\n5KtknANcBjwhaXW27Xrgw5IWka4yeQ64ojlNNLN24EqGmdUrz9Ulv6ByxeOHY98cM2tXecdkeOCn\nmZV4xk8zy8WVDDOrl0OGmeWSZ0xGaZ6MUiUjYnzbaGbtxSHDzHLJW8no7k73N+nudpeJ2WTnkGFm\nudQzJgM8LsPMHDLMLKfhlYyRuku6utK6x2WYmUOGmeVSHjK6uyvfm6S/f7Da4ZBhZg4ZZpZLeciY\nNs0hw8xqc8gws1zKx2SMNKhz377B28F7TIaZOWSYWS6uZJhZvfLcu2ShpIckPSnpV5I+nW0/WtJK\nSU9LesC3ejeb2IaPyahVyXDIMLM8lYx+4JqIeAtwNvCXkk4BrgVWRsRJwIPZczOboHbvhpkz03qe\nSsasWbBr1/i1z8zaT82QERG9EbEmW98BPAUcB1wMLM8OWw5c0qxGmlnrlYeMPJWMww6DnTvHr31m\n1n7qGpMh6Q3AYuBhYF5E9GW7+oB5Y9oyM2sru3en6gRUv4S1FDJcyTCzPLd6B0DSbOAu4OqI2C7p\n4L6ICEkV71LQ09NzcL1QKFAoFBptq5m10K5dQ7tLRqpklLpLXMkw63zFYpFisdjw63OFDEndpIBx\ne0SsyDb3SZofEb2SFgCbK722PGSYWeca3l1Sq5LhkGHW+YYXB5YtW1bX6/NcXSLgFmBdRHy9bNc9\nwNJsfSmwYvhrzWzi2LVrsLskz8BPhwwzy1PJOAe4DHhC0ups23XAl4E7JV0ObAA+2JQWmllbqHfg\n56xZDhlmk13NkBERv2Dkisf5Y9scM2tXtS5hjTi0u2TLlvFto5m1F8/4aWa5lHeXVKpkHDgA0uBd\nWN1dYmYOGWaWS61KRvl4DHDIMDOHDDPLqdaYjPLxGOAxGWbmkGFmOe3cmaoTkL+S4cm4zCY3hwwz\nqykCtm+H2bPT8zyVDHeXmJlDhpnVtHdvGtBZqlRUmvGz/MoScMgwM4cMM8th+3Y4/PDB55VCRvmt\n4MFjMszMIcPMctixY2jImDED9uwZesyePWl7icdkmJlDhpnVVD4eA/KHDFcyzCY3hwwzq6lSJWPv\n3qHH7N3rkGFmQ+W5QdqtkvokrS3b1iNpk6TV2bKkuc00s1YaXsmYPj1VLiIGtw2vZJSCyIED49dO\nM2sveSoZtwHDQ0QAN0bE4mz50dg3zczaxfCBn11dMHXq0LkyhoeMKVM8+NNssqsZMiLi58CrFXZp\n7JtjZu1o61aYM2fotuHjMoaHDICjjkqvNbPJaTRjMq6S9LikWyTNqX24mXWqrVtTYCiXJ2QcfTS8\n/HLz22dm7anmrd5HcBPwhWz9i8BXgcsrHdjT03NwvVAoUCgUGvxIM2uVV19trJLxutfBK680v31m\n1hzFYpFisdjw6xsKGRGxubQu6Wbg3pGOLQ8ZZtaZtm6F179+6LbS4M+SPXuGTsYFrmSYdbrhxYFl\ny5bV9fqGukskLSh7+n5g7UjHmlnncyXDzBpRs5Ih6Q7gPcBcSRuBzwMFSYtIV5k8B1zR1FaaWUs1\nOvDTlQyzya1myIiID1fYfGsT2mJmberll1NgKDd8Qq7hk3FBqmS8+GLz22dm7ckzfppZTVu2wDHH\nDN02Ywbs3j34fPduj8kws6EcMsyspkohY/hEWzt3pqnEy3lMhtnk5pBhZlXt2ZNu637EEUO3z559\naMgon3ocUiXDIcNs8nLIMLOqXnoJ5s4FDZvjd/gN0EaqZLi7xGzycsgws6p6e+HYYw/dnidkHHMM\n9PU1t31m1r4cMsysqk2bYOHCQ7fnCRlz58KuXWkxs8nHIcPMqtq0CY4//tDthx0GO3YMPt+x49CQ\nIcFxx8FvftPcNppZe3LIMLOqNm5svJIBKaBs2tS89plZ+3LIMLOqqlUy8oSMhQsdMswmq5ohQ9Kt\nkvokrS3bdrSklZKelvSAb/VuNnGNFDLyXMIK6bUbNzavfWbWvvJUMm4Dlgzbdi2wMiJOAh7MnpvZ\nBDRSd8ns2bB9e1qPcHeJmR2qZsiIiJ8Drw7bfDGwPFtfDlwyxu0yszYwMAC//W0avDncnDmwbVta\n37EjTTPe3X3oca5kmE1ejY7JmBcRpavf+4B5Y9QeM2sjL7yQLkOdOfPQfXPmpLuzQuW7tJa86U3w\nzDPNa6OZta9RD/yMiCDd8t3MJpj16+GUUyrvO/LIwUrG1q3peSUnnggbNkB/f1OaaGZtrOat3kfQ\nJ2l+RPRKWgBsHunAnp6eg+uFQoFCodDgR5rZeFu/Hk4+ufK+8krGtm0jVzKmT09jOn7965Hfy8za\nU7FYpFgsNvz6RkPGPcBS4IbsccVIB5aHDDPrLOvXw+mnV9532GGwd2+qUFTrLoFUDakWWMysPQ0v\nDixbtqyu1+e5hPUO4P8Cb5a0UdIngC8DF0h6Gjgve25mE8xTT40cDKTBaka17hJI7/HUU81po5m1\nr5qVjIj48Ai7zh/jtphZGxkYgCeegNNOG/mY0l1WS3dqHcnpp8N99419G82svXnGTzOr6NlnU3Vi\nXpVrx+bNS3dZ7eurftyZZ8KqVWPfRjNrbw4ZZlbRI4+kcFBN3pBx4onwyiuwZcvYttHM2ptDhplV\ntGpVvpDR21s7ZEyZAmec4WqG2WTjkGFmFf3sZ3DOOdWPmT8fXnwxBY3586sfe8458G//NnbtM7P2\n55BhZofo7YXnn0/Vh2p+7/fS2I1nnknr1VxwATzwwNi10czan0OGmR3iJz+Bc8+FqTWuP3vzm1PF\nY/p0OPro6seeeSY891zqWjGzycEhw8wO8a//ChdeWPu4E09MgzkXLap9bHd3qmbce+/o22dmncEh\nw8yG2LEjdWv8yZ/UPvbww+Fv/xauuirfe3/oQ/DP/zy69plZ51C6v1mT3lyKZr6/mY2922+HO+6A\n++8f+/fetSvdNv7JJ+H1rx/79zez5pJERCjv8a5kmNkQ//iP8Od/3pz3njULPvIRuOmm5ry/mbWX\nUVUyJG0AXgMOAP0Rceaw/a5kmHWQRx9N3ST/+Z+1B3026j/+A979bnjhBZgxozmfYWbNMd6VjAAK\nEbF4eMAws87zd38H117bvIAB6YqUM86AW25p3meYWXsYbSXjOeAdEfHyCPtdyTDrEL/4BVx2Wao0\nTJ/e3M9aswaWLEmfVe3urWbWXlpRyfiJpEclfXKU72VmLbJvH3zqU/ClLzU/YEC65PXCC+ELX2j+\nZ5lZ64y2KHpORLwo6RhgpaT1EfHz8gN6enoOrhcKBQqFwig/0szG2t//Pfzu78Kll47fZ95wQwob\nF18M73nP+H2umeVXLBYpFosNv37MLmGV9HlgR0R8tWybu0vM2twPfwh/9mfprqvjfVnp/ffDX/wF\nPPwwLFgwvp9tZvUbt+4SSbMkHZ6tHwb8EbC20fczs/H32GOwdGmaIKsV81ZceCFccQX88R/Da6+N\n/+ebWXONZkzGPODnktYADwM/iAjf/sisQzz2WPoh/93vwrve1bp2XH89nH02nH8+vPJK69phZmPP\nM36aTUJ33ZW6Kb77Xbjkkla3BiLgc5+DH/wA7r4bTj651S0ys0rq7S5p4tXwZtZuduyA666De+6B\nH/0I3v72VrcokeArX0lzaLz73Wn94x9P282sc3lacbNJIAL+5V/gtNNS0Fizpn0CRrnLL4eVK+Ef\n/gHOOy/d48TMOpdDhtkENjAA992XxjwsWwbf+Q7cdhscdVSrWzayRYvS1Sbvf38KGh/5iMOGWafy\nmAyzCai3N91J9dvfTjNq/vVfwwc/CFM67NeK7dvhm99My8knw5VXpnk1fM8Ts9aod0yGQ4bZBLFh\nQxpnceedsHp1uiz0yitTFaPTxzbs2wcrVsA//dPgVTF/+qfw3vd6WnKz8eSQYTYJRMBzz8GqVeme\nIw88ANu2wQUXwAc+kO4LMnNmq1vZHH19aXzJ3XfDL3+Zxpm8973pMtwzzoDXva7VLTSbuBwyzCaY\nXbvSjcSeegrWrUuDNletgu5uOOsseOc7U7h461s7rztktPbsSUHjJz9Jj48+Cscck8LG4sVwyilw\n6qlwwgnQ1dXq1pp1PocMsw5z4EAaQ/H882l54YX0uGEDrF8PL74IJ544+APz9NNTuDjuuFa3vP0M\nDKRAtmoVPPHEYDDbvBne9CY46aR0j5bhy5w5nd+lZDYeHDLMWmxgIF0munUrvPoqbNmSfsj19Q19\nLK339qarPSr98Hvzm+GNb4SpntFmVHbuTOHj2WcHw1wpyD3/fOp+mjcvLfPnD66Xlrlz09/RUUel\nQDJrlkOJTU4OGWZ1GhiA3bvTD6KdO1P3RK31115LIaLS8tpr6YfQnDlpUOKxx6Zl3rzKjwsW+GqJ\nVopIf2elwNfXd+jy8sspML76avo7PnAg/f2WgkcpfBx+OMyePXQ57LBDt5W2z5yZ/u67ux1arDOM\na8iQtAT4OtAF3BwRNwzb75AxRorFIoVCodXNqCoi/cDevx/6+9PjSEu1/dX27dsHe/cOPuZdhh//\n2mtFoMDOnalff8aM9KU/a1Z6HL4+/PkRR6QfKuXLkUcOPk62ykMn/PscS3v2DFaqysPHjh0jLzt3\nDn2+fXt6nz170v+bGTMGl4giRx1VGLKt0tLdPbhMmzb0eaPbu7oGl6lThz4faWnngDTZ/m0227hN\nKy6pC/gWcD7wG+ARSfdExFONvmcr3HwzfPKTrW5FHkWg0OI2tLfp0ysv06alxyOOGNz2618XWby4\ncPC3ybwDJvfvT1dxbNsGGzc298/TSX75yyLvfGeh1c1oO6UwMHdu9eP270/htxQ6Hn+8yHHHFQ4+\nLw8kpWX37vQ6q6VIu393nn9+mul2IhrN71tnAs9GxAYASf8HeB/QUSFj/nz4nd9J69JgIh++Xmt/\nnvXRvMfGjUPbOWVK+g1iypT8Sz3H1/vepd94OuXqhn374B3vaHUrJo516wb/fdrobduWLkPuBKUK\n5oEDg8v+/UOfV1tKrx/NY7V9GzbAwoX1vWf5n63SY95teY8/9dRR/RW0tdGEjOOA8t/lNgFnja45\n4++ii9LS7np60mJjY9s2uOaaVrdi4vD5HFs+n2PH352t1fCYDEkfAJZExCez55cBZ0XEVWXHeECG\nmZnZBDJet3r/DbCw7PlCUjWjoYaYmZnZxDKaHvRHgRMlvUHSNOBDwD1j0ywzMzPrdA1XMiJiv6S/\nAn5MuoT1lk67ssTMzMyap6mTcZmZmdnk1ZQLDiVdJekpSb+SdEPZ9uskPSNpvaQ/asZnT1SS/quk\nAUlHl23z+ayTpP+e/dt8XNLdko4s2+fzWSdJS7Lz9Yykz7W6PZ1G0kJJD0l6Mvu+/HS2/WhJKyU9\nLekBSXNa3dZOIalL0mpJ92bPfS4bJGmOpO9n35nrJJ1V7/kc85Ah6VzgYuCtEXEa8D+y7aeSxm2c\nCiwBvi2pQ2ZVaC1JC4ELgOfLtvl8NuYB4C0R8fvA08B14PPZiLIJ+ZaQztuHJZ3S2lZ1nH7gmoh4\nC3A28JfZObwWWBkRJwEPZs8tn6uBdUCpTO9z2bhvAPdHxCnAW4H11Hk+m/EleiXwpYjoB4iILdn2\n9wF3RER/NoHXs6QJvay2G4HPDtvm89mAiFgZEaXpdh4Gjs/WfT7rd3BCvuz/e2lCPsspInojYk22\nvoM0meFxpF/UlmeHLQcuaU0LO4uk44ELgZuB0tWNPpcNyKq8746IWyGNw4yIbdR5PpsRMk4E/lDS\nv0sqSirNq/h6hl7iuon0n8mqkPQ+YFNEPDFsl8/n6P0X4P5s3eezfpUm5PM5a5CkNwCLSeF3XkT0\nZbv6gHktalan+RrwGaBs3k6fywadAGyRdJuk/yfpu5IOo87z2dDVJZJWAvMr7Pqb7D2PioizJZ0B\n3Am8cYS38qhTap7P64Dy8QHV5h7x+aTq+bw+Ikr9tH8D7IuI/13lrXw+q/P5GSOSZgN3AVdHxHaV\n3XEsIsITG9Ym6SJgc0SsllSodIzPZV2mAm8D/ioiHpH0dYZ1jeQ5nw2FjIi4YKR9kq4E7s6OeyQb\nrDiXQyfvOj7bNumNdD4lnUZKk49nXzrHA49JOgufzxFV+/cJIOnjpJLqe8s2+3zWr+aEfFabpG5S\nwLg9IlZkm/skzY+IXkkLgM2ta2HH+APgYkkXAjOAIyTdjs9lozaRquiPZM+/T/qlt7ee89mM7pIV\nwHkAkk4CpkXES6SJui6VNE3SCaRulVVN+PwJIyJ+FRHzIuKEiDiB9Jf+tqxU5fPZAElLSOXU90XE\nnrJdPp/184R8o6T028MtwLqI+HrZrnuApdn6UtL3qlUREddHxMLsu/JS4KcR8VF8LhsSEb3Axuzn\nOKQ7rj8J3Esd53M004qP5FbgVklrgX3Ax7IGr5N0J2nU737gU+FJOup18Hz5fDbsm8A0YGVWHfpl\nRHzK57N+npBvTJwDXAY8IWl1tu064MvAnZIuBzYAH2xN8zpa6f+vz2XjrgL+V/ZLxK+BT5D+r+c+\nn56My8zMzJrC8wCYmZlZUzhkmJmZWVM4ZJiZmVlTOGSYmZlZUzhkmJmZWVM4ZJiZmVlTOGSYmZlZ\nU/x/8QPGQJ9lUrYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xad52978>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# fourier transform\n",
    "N = len(t)\n",
    "dt = t[1]-t[0]\n",
    "\n",
    "# calculate the fast fourier transform\n",
    "# y2 is the solution to the under-damped oscillator from the previous section\n",
    "F = fft(y2[:,0]) \n",
    "\n",
    "# calculate the frequencies for the components in F\n",
    "w = fftfreq(N, dt)\n",
    "\n",
    "fig, ax = subplots(figsize=(9,3))\n",
    "ax.plot(w, abs(F));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear Algebra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A = array([[1,2,3], [4,5,6], [7,8,9]])\n",
    "b = array([1,2,3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.33333333,  0.66666667,  0.        ])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# solve a system of linear equations\n",
    "x = solve(A, b)\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.34211023,  0.22314378, -0.11368553])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# eigenvalues and eigenvectors\n",
    "A = rand(3,3)\n",
    "B = rand(3,3)\n",
    "\n",
    "evals, evecs = eig(A)\n",
    "\n",
    "evals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.24560175, -0.93054452,  0.46758953],\n",
       "       [-0.84521949, -0.00332278, -0.56599975],\n",
       "       [-0.4746407 ,  0.36616371,  0.67897299]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "evecs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[-0.23411483,  0.96016594,  0.15255034],\n",
       "        [-0.84815445, -0.12501447, -0.51478677],\n",
       "        [-0.47520972, -0.24990547,  0.84363676]]),\n",
       " array([ 1.3457633 ,  0.23442523,  0.1079208 ]),\n",
       " array([[-0.21657397, -0.81553512, -0.53665463],\n",
       "        [ 0.7658064 ,  0.19902251, -0.61149865],\n",
       "        [-0.60550497,  0.54340824, -0.58143892]]))"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "svd(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy import optimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEACAYAAACnJV25AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHxdJREFUeJzt3XmYVNWZx/HvyyIIKMgiIIsQBSKKKxIRI+2SiIlCYozG\nzIxGjfGJjMk4WUazIJk8Gk2iMTGT5YkGTVSMS1zQqCjSEREhJKBCQwABpYk2CDayr+/8cartsu2G\nqu7qOrdu/z7Pcx9quVX3paF/dercc841d0dERNKjVewCRESksBTsIiIpo2AXEUkZBbuISMoo2EVE\nUkbBLiKSMjkFu5n1M7PpZrbQzBaY2dcyj3c1s2fNbImZTTWzLlmvuc7MlprZYjP7ZHP9BURE5IMs\nl3HsZtYL6OXu882sE/B34DPApcA77v5jM/sf4CB3v9bMhgL3AScCfYDngMHuvqe5/iIiIhLk1GJ3\n97fdfX7m9iZgESGwxwJ3Z3a7mxD2AOOAye6+091XAsuAEQWsW0REGpB3H7uZDQCOA2YDPd29KvNU\nFdAzc/sQoDLrZZWEDwIREWlmeQV7phvmYeDr7r4x+zkPfTp769fR2gUiIkXQJtcdzawtIdT/6O6P\nZh6uMrNe7v62mfUG1mQeXw30y3p538xj2e+noBcRaQR3t709n+uoGAPuBCrc/baspx4HLsncvgR4\nNOvxL5jZfmY2EBgEzKmnuMRv119/ffQaVKfqVJ2qsWbLRa4t9lHAvwOvmtm8zGPXATcBD5jZ5cBK\n4IJMYFeY2QNABbALuMpzrUhERJokp2B39xdpuHV/ZgOvuRG4sZF1iYhII2nm6T6UlZXFLiEnqrOw\nVGdhlUKdpVBjrnKaoNQsBzZT74yISJ7MDC/EyVMRESkdCnYRkRKxJ8dFWRTsIiIl4sUXc9tPwS4i\nUiJmzcptPwW7iEiJyDXYNSpGRKQEuEOvXrBmjUbFiIikwooV0CbHtQIU7CIiJWDWLBg5Mrd9Fewi\nIiVAwS4ikjL5BLtOnoqIJNzmzXDwwbBuHey/v06eioiUvLlzYdgwaN8+t/0V7CIiCZdPNwwo2EVE\nEu/ll/MLdvWxi4gkWM3EpLlzoV8/LdsrIlLyVqyAtm1DqOcqarCrwS4isnf59q9D5GBfsSLm0UVE\nkm/WLDjppPxeEzXYX3op5tFFRJJvxgw45ZT8XqNgFxFJqOpqWL4cjj8+v9cp2EVEEmrmTBgxIpw8\nzUfUYF+6FN57L2YFIiLJNWMGfPzj+b8uarAffzzMmROzAhGR5CrJYD/5ZHXHiIjUZ+tWmD8//xEx\nkIBgnzkzZgUiIsk0Zw4cdRR07Jj/a6MG+8iRYQ2E3btjViEikjyN7YaByMF+8MFhq6iIWYWISPKU\nbLCD+tlFROratSv0ZuQ7MalG9GAfNUrBLiKS7ZVXoG9f6Natca+PHuxqsYuIfFBTumEgAcE+dCis\nXQtr1sSuREQkGV54ocSDvVWrME5TrXYRkbCc+YsvlniwQzhBMGNG7CpEROKrqAhj1/v3b/x7JCLY\nR48OXz1ERFq66dPh9NOb9h45BbuZ/d7MqszstazHJppZpZnNy2xnZz13nZktNbPFZvbJfb3/iSfC\nokWwcWPj/hIiImnx/PNFCnZgEjCmzmMO3Orux2W2pwDMbChwITA085pfmdlej9O+PZxwgpYXEJGW\nbc8e+Otf4bTTmvY+OQW7u88A3q3nqfqulD0OmOzuO919JbAMGLGvY6g7RkRauldegR494JBDmvY+\nTe1jv9rMXjGzO82sS+axQ4DKrH0qgT77eqNTT1Wwi0jLNn1601vrAG2a8NpfA/+buf1D4Bbg8gb2\n9foenDhx4vu3P/axMubPL2PrVth//yZUJSJSop5/Hi655IOPlZeXU15entf7mHu9mfvhHc0GAFPc\nfdjenjOzawHc/abMc08D17v77Dqv8brHHjkSbryxMJ9YIiKlZNeusITAsmWhO6YhZoa719cN/r5G\nd8WYWe+su58FakbMPA58wcz2M7OBwCAgp+skqTtGRFqqv/8dDj1076Geq5y6YsxsMjAa6G5mq4Dr\ngTIzO5bQzbICuBLA3SvM7AGgAtgFXPWhpnkDTj0Vbr01/7+EiEipK1T/OuTRFVNo9XXFbNgQVjRb\ntw722y9KWSIiUXzykzB+PIwbt/f9mrUrpjl07gyDBsHcubErEREpnh07wvrro0cX5v0SFeygfnYR\naXnmzIHBg6FLl33vm4vEBfvo0WHmlYhIS/Hcc01fRiBbIoN95szw1UREpCV45hk466zCvV/igr1r\n19DPPienAZIiIqXt3Xdh4cLGX9+0PokLdoAzzoBp02JXISLS/KZNC6Herl3h3lPBLiIS0dSpYahj\nISVqHHuNzZuhZ0+oqgpXEhERSSN3GDAAnn4ajjgit9eU3Dj2Gh07hvXZdbk8EUmzJUvCGuwf/Whh\n3zeRwQ7qjhGR9KsZDWN7bX/nT8EuIhJJc/SvQ0L72AF27oTu3WH58rCUpYhImmzfHlZyXLkyDPPO\nVcn2sQO0bRuGAE2fHrsSEZHCe+mlcMI0n1DPVWKDHdQdIyLp1VzdMKBgFxGJ4umnC7uMQLZEB/uw\nYWG67Ztvxq5ERKRwKitDrp10UvO8f6KDvVWr8FXlmWdiVyIiUjhPPglnnw1tcrqGXf4SHewAY8aE\nrywiImnxxBPw6U833/sndrhjjaoqGDIE1q4NI2VERErZli3Qq1f+wxxrlPRwxxo9e8Jhh8GsWbEr\nERFpuunT4bjjmmeYY43EBzuEvih1x4hIGjz5JJxzTvMeoySCXf3sIpIG7qF/XcFOGBK0YgW8/Xbs\nSkREGu+118K5wkKv5lhXSQR7mzZhstLUqbErERFpvJrWeqFXc6yrJIId1B0jIqWvGN0wUALDHWus\nWhXOJFdVQevWzViYiEgzWLsWDj8c1qxp2vVNUzHcsUa/fmHs59y5sSsREcnf44+HtWEKedHqhpRM\nsEPojnnqqdhViIjk789/hvPOK86xSirYzz0XpkyJXYWISH42bAjXcP7Up4pzvJIK9lGjwjTcysrY\nlYiI5O4vf4FTT4UDDyzO8Uoq2Nu0CZ94arWLSCkpZjcMlFiwA4wdG05CiIiUgi1bwhycsWOLd8yS\nC/azzoIXX4SNG2NXIiKyb1OnwvDh0L178Y5ZcsF+4IFw8smahSoipaHY3TBQgsEO6o4RkdKwY0eY\nbfqZzxT3uDkFu5n93syqzOy1rMe6mtmzZrbEzKaaWZes564zs6VmttjMCn4d7nPPDUtf7tpV6HcW\nESmc8vJwoaA+fYp73Fxb7JOAMXUeuxZ41t0HA9My9zGzocCFwNDMa35lZgX9ZtC/f5iJqotviEiS\nPfRQ8bthIMdgd/cZwLt1Hh4L3J25fTdQ82VjHDDZ3Xe6+0pgGTCi6aXWObi6Y0QkwXbsCP3rF15Y\n/GM3pSXd092rMrergJ6Z24cA2VOIKoGCfxEZOxYeeywsXC8ikjRTp4Z11/v3L/6xC9JFklmmcW8R\nW/D4Pf542L4dFiwo9DuLiDTd5Mlw0UVxjt2mCa+tMrNe7v62mfUG1mQeXw30y9qvb+axD5k4ceL7\nt8vKyigrK8v54GZw/vmhD2vYsDwrFxFpRlu2hAEeP/tZ09+rvLyc8vLyvF6T83rsZjYAmOLuwzL3\nfwysc/ebzexaoIu7X5s5eXofoV+9D/AccHjdxdfzXY+9Pi+/DJddBhUVTXobEZGC+tOf4Pe/h2ee\nKfx7F2w9djObDLwEDDGzVWZ2KXAT8AkzWwKcnrmPu1cADwAVwFPAVU1O8AaMGBFmoCrYRSRJ7r8f\nvvCFeMcvmSsoNeSaa+Cgg2DChAIUJSLSRNXVcOih8MYb0KXLvvfPV6quoNSQ88+HBx+MXYWISPDI\nI3D66c0T6rkq+WAfORLWr4fFi2NXIiISdzRMjZIP9lat4HOfC6NjRERiWr0a/vY3OOecuHWUfLBD\n7bBHEZGY7rkn5FGHDnHrSEWwjxoFVVWwdGnsSkSkpXKHu+6CL30pdiUpCfbWreHznw99WyIiMcye\nDbt3h+tFxJaKYAf4t3+De+/V2jEiEkdNa932OhCxOFIT7CNGwJ49MHdu7EpEpKXZujUMu7744tiV\nBKkJdrPaVruISDE9+mi4rmnfvrErCVIT7BCC/f77dWUlESmuu+6CSy+NXUWtVAX7oEFhKu+0abEr\nEZGWorIydAGPGxe7klqpCnYIrfZ77oldhYi0FHfeGa6StP/+sSupVfKLgNVVVRUuHrt6NXTsWPC3\nFxF5365dMGAAPPVU8a4L0SIWAaurZ8+wfsxjj8WuRETSbsqUEOxJu9hP6oId4D/+A/7wh9hViEja\n/frX8NWvxq7iw1LXFQNhTGnfvjBvXpwLyYpI+i1dGpYzWbUK2rUr3nFbZFcMhJMYF10EkybFrkRE\n0uq3vw0zTYsZ6rlKZYsdYP78MPxo+fKwloyISKFs3Rp6A15+GQ47rLjHbrEtdoBjj4Xu3TWmXUQK\n78EH4YQTih/quUptsAN8+ctwxx2xqxCRNHGHX/4SrroqdiUNS21XDISLyg4YEE5y9OjRrIcSkRbi\nxRfD8gGLF8fp5m3RXTEQLiY7dqxmoopI4dxyC1xzTbLP3aW6xQ7wwgthnOmCBclYJ1lESteyZWEC\n5MqV8Wa2t/gWO8DHPx7WaX/hhdiViEip+/nP4StfSf5yJalvsUM40fHXv4Yz2SIijbF+PRx+OCxc\nCL17x6sjlxZ7iwj2jRvDcr6vvAL9+hXlkCKSMjfdFE6Y3nVX3DoU7Fm+9jU44AC44YaiHVJEUmL7\ndvjIR8IqjkcfHbcW9bFnGT8+jGnfti12JSJSau6+OwR67FDPVYsJ9iFDwmzUBx6IXYmIlJKdO+FH\nP4IJE2JXkrsWE+wAV18Nt98eZo6JiOTinnvC0gEjR8auJHctKtjPPhvWrYPZs2NXIiKlYNeucF7u\n+9+PXUl+WlSwt24dTqLeckvsSkSkFNx/P/TpA6NHx64kPy1mVEyNTZtg4EB46SUYNKjohxeRErF7\nNxx5ZJgHc+aZsauppVEx9ejUKSwx8NOfxq5ERJLswQfhoIPgjDNiV5K/FtdiB1i7NoySqaiAXr2i\nlCAiCbZzJwwdGq5pmqTWOqjF3qAePeCLXwzrPoiI1HXHHaHLNmmhnqsmt9jNbCXwHrAb2OnuI8ys\nK/An4FBgJXCBu1fXeV20FjvAihUwfHi4dF7nztHKEJGE2bQJBg+GJ56A44+PXc2HFavF7kCZux/n\n7iMyj10LPOvug4FpmfuJMnAgjBkTLkgrIlLjttvCKJgkhnquCtFiXwEMd/d1WY8tBka7e5WZ9QLK\n3f2jdV4XtcUO8OqrcNZZ8Prr0KFD1FJEJAHWroUjjghzXZJ6PdNittifM7O5ZnZF5rGe7l6VuV0F\n9CzAcQru6KNh1Cj41a9iVyIiSXDDDXDRRckN9VwVosXe293fMrMewLPA1cDj7n5Q1j7r3b1rnddF\nb7FDWFv59NPDlVEOOCB2NSISy+LF4cI8CxZAz0Q2RYNcWuxtmnoQd38r8+daM3sEGAFUmVkvd3/b\nzHoDa+p77cSJE9+/XVZWRllZWVPLyduRR4ZxqrffDt/5TtEPLyIJ4B7Wkvre95IX6uXl5ZSXl+f1\nmia12M2sA9Da3TeaWUdgKvAD4ExgnbvfbGbXAl3c/do6r01Eix3gn/+EU04JrXaNkBFpef78Z7j+\nepg3D9o0ubnbvJr9QhtmNhB4JHO3DXCvu/8oM9zxAaA/CR3uWNeXvgQDBkDWlwgRaQG2bAmTke66\nCyJ0GuRNV1DKw/LlMGJEaL136xa7GhEplgkTYMmSsOBXKVCw52n8eGjbNoxjFZH0W7YMTjoJ5s+H\nvn1jV5MbBXue1q4NX8lmzgwzz0QkvfbsCSPixo6F//7v2NXkTmvF5KlHD/j2t+Fb34pdiYg0t9/8\nJlyk+utfj11J4anFXse2baHVfscd4dNcRNJn5cqwVtSMGWGmaSlRi70R2reHm2+Gb3wjLLQvIuni\nDl/5SvgdL7VQz5WCvR7nnw8dO4bhTyKSLnfeGa59/M1vxq6k+agrpgH/+Ee4+PXChdC9e+xqRKQQ\naiYjPv88DBsWu5rG0aiYJrrmGqiuhkmTYlciIk21fTuMHAlXXBEuj1mqFOxNtHFjWEvmD38ojRlp\nItKw//ovePNNePhhsL3GYrIVZRGwNDvggLA42JVXhrXb27WLXZGINMYTT8Ajj4S1YEo51HOlk6f7\nMG5cGP54002xKxGRxli5Er78ZbjnHujadZ+7p4K6YnKwalW4TNZzz8Exx8SuRkRytXkznHwyXHpp\n6IpJA/WxF9Ddd8Mtt8CcOWGsu4gkmztccEEYujxpUnq6YDRBqYAuvhgGDYLvfz92JSKSixtuCN+2\nf/Ob9IR6rtRiz8M774SumHvv1SgZkSR7+OGwBszf/ga9e8euprDUYi+w7t3hd78LF+XYsCF2NSJS\nn+nTwzj1KVPSF+q5Uou9Ef7zP+Gtt+Chh1reVzyRJPvHP2DMGHjggfR+q1aLvZnccgtUVsJPfxq7\nEhGpsXQpnHMO/Pa36Q31XKnF3khvvhkupXf//fpPJBLbsmVw5pnwve+FMetpphZ7M+rfP0x4+OIX\nYfXq2NWItFyLFoXG1Xe+k/5Qz5WCvQnOPDP0t593XrjSuYgU16uvwhlnwI03hjXWJVBXTBO5wyWX\nwHvvhSFWrVvHrkikZZg5Ez73Ofj5z+HCC2NXUzzqiikCs3AZvffeC1OWU/BZJZJ4994Ln/1suBhO\nSwr1XKnFXiDV1WEB/8suK60rnouUEnf4wQ/CEh9TpsBRR8WuqPi0bG8RdekCf/kLjBoVbl92WeyK\nRNKlujqcHK2shJdfhp49Y1eUXOqKKaD+/WHaNJgwQVddEimkOXPCCqu9e0N5uUJ9X9RiL7DBg0O4\nn346tGoVTqyKSOPs3g0/+xn8+MdhMa/zzotdUWlQsDeDIUNCuJ9xBuzaBZdfHrsikdKzYEHoetlv\nP5g9GwYOjF1R6VBXTDP56EfDYkQ33AATJ2q0jEiutm0L3ZmnnRbOVZWXK9TzpWBvRoMHw6xZ4XqL\nl18OO3fGrkgkufbsgT/+MXzjXbgQ5s8Pk45aKaXypuGORbBpUxhru3MnTJ4M3brFrkgkOdzh6afD\nkgDt28NPfhKGDkv9NEEpITp1gsceg6OPhhNOCP2FIi3drl1w331w3HHw7W+HYH/pJYV6IajFXmSP\nPAJXXhn6EMeP13ru0vJUVoYZo3fcAYceGkL9U5/S70KudDHrhHr9dfj856FXr7B2dL9+sSsSaV7V\n1eFc0733hm+sF14YzjsNHx67stKjYE+wHTvg5pvhF78II2euuEItlubmDuvXQ1VVuH5tzVZdHbYN\nG8KaP5s3h23LFti+Pfxbbd8eug727Ambe/j3atUqbG3bhmF57dqFbf/9oUOHsB1wQO3WuXOYmdyl\nCxx0EHTtGrZu3cLr0sI9LKf7/PPw+ONhpmhZGVxwQRiL3qFD7ApLl4K9BCxYEIZ0tW8frsx04omx\nKypN7rBmDaxcCW+8Ef6srKzd/vWv8HynTnDwwdCjR7iGbbduIWA7dw7bgQdCx45h69Ah/Lu0axdC\nu23b2iCvOaZ7mESza1f4ANixA7Zurd02b4aNG2u3996r/SBZvz5s69aFP9u1CzV17x7qq6nz4IPD\nTMvsP3v0CDUlxbvvwrx5YZszJwxR7NgRRo+GT386XK6uU6fYVaZD1GA3szHAbUBr4A53v7nO8wr2\njN27Q5/jhAmhVXPjjaHvUT7IPVxrdsmS2m3ZstC1tXx5aCUPHBh+doceGpZ46NMH+vaFQw4JoZjU\nVrF7CP533oG1a2u3NWvCN4y1a8OfNbfXrg1B2aPHBz8EunWr3bp2rf1m0Llz+MbQqVN+Hwju4dvK\nhg3hA+itt8K2enX4uS9bFi5JV10NxxwTpv0PHw6nngoDBjTbj6tFixbsZtYa+CdwJrAa+Btwkbsv\nytpHwV7Hpk3hOqq33x6+rn7zm2FMb0uzezesWAEVFWFbtChsixeHFvSQIWGOwKBBYTvsMPjIR0Jr\nu6XYsyeE6Zo1H/wgWLcufDisWxda0TXfDqqrw/+vjRtDF1L79rXfRtq2re1WgjAst6YLqmb/zp3D\nB0Xv3rXbYYeFn//hh4cPUV2LoDhiBvtI4Hp3H5O5fy2Au9+UtY+CvQHvvAP/939hO+UUuPrq8JU2\nbRM13GHVqtAdlb0tXhxaoEOHwpFHwhFHhG3IkBAu0njuIbC3bQvhvW1bCPKabiX30KKv2Tp1Ch8A\nkhwxg/184Cx3vyJz/9+Bj7n71Vn7KNj3YfPmsErk734X+mYvuQQuvji0TkuJO7z9dphNuHBhbYBX\nVITgOOqoEOBHHRW2oUPVHyvSkJjrseeU2BMnTnz/dllZGWVlZc1UTmnq2DFcU3X8+HBSatIkOOmk\n0Jo999ywnXhi+CqdBHv2hBOVNV0nixaF8F64MHzbGDoUhg0L/bAXXxzCXC1wkb0rLy+nvLw8r9c0\nV4v9JGBiVlfMdcCe7BOoarE3zu7dYdTBlCnw5JPhpOHw4eECHyeeGMJy4MDm6+/cujWE9xtv1J60\nfP312hOZXbqEBdBquk9qulMOPljDOUUKIWZXTBvCydMzgH8Bc9DJ02ZRXR0WGps5M7TqKyrCyIlB\ng8LEp759w8iQbt1qh/R17BiCv02b0JKuGae9fXs4wVYzprtmzHfNVlkZHu/bN5wsO+yw2hOXgweH\nk2gHHBD7JyKSbrGHO55N7XDHO939R3WeV7A3k02bwhC07HHc69eHUN6wIfTd14y93rPngxNrOnUK\nre7OncMwuZ49a7e+fUPLO20ncUVKiSYoiYikjFZ3FBFpgRTsIiIpo2AXEUkZBbuISMoo2EVEUkbB\nLiKSMgp2EZGUUbCLiKSMgl1EJGUU7CIiKaNgFxFJGQW7iEjKKNhFRFJGwS4ikjIKdhGRlFGwi4ik\njIJdRCRlFOwiIimjYBcRSRkFu4hIyijYRURSRsEuIpIyCnYRkZRRsIuIpIyCXUQkZRTsIiIpo2AX\nEUkZBbuISMoo2EVEUkbBLiKSMgp2EZGUUbCLiKSMgl1EJGUU7CIiKaNgFxFJmUYHu5lNNLNKM5uX\n2c7Oeu46M1tqZovN7JOFKVVERHLRlBa7A7e6+3GZ7SkAMxsKXAgMBcYAvzKzkv1mUF5eHruEnKjO\nwlKdhVUKdZZCjblqauBaPY+NAya7+053XwksA0Y08TjRlMo/tuosLNVZWKVQZynUmKumBvvVZvaK\nmd1pZl0yjx0CVGbtUwn0aeJxREQkR3sNdjN71sxeq2cbC/waGAgcC7wF3LKXt/LClSwiIntj7k3P\nXDMbAExx92Fmdi2Au9+Uee5p4Hp3n13nNQp7EZFGcPf6usHf16axb2xmvd39rczdzwKvZW4/Dtxn\nZrcSumAGAXPyLUxERBqn0cEO3GxmxxK6WVYAVwK4e4WZPQBUALuAq7wQXwtERCQnBemKERGR5Ig6\nvryeSU5jYtazL2b2DTPbY2ZdY9dSHzP7YWaU0nwzm2Zm/WLXVB8z+4mZLcrU+mcz6xy7pvqY2efN\nbKGZ7Taz42PXk83MxmQmAC41s/+JXU99zOz3ZlZlZq/te+94zKyfmU3P/FsvMLOvxa6pPmbW3sxm\nZ36/K8zsRw3tG3viUN1JTk9HrqdBmZD8BPBG7Fr24sfufoy7Hws8Clwfu6AGTAWOdPdjgCXAdZHr\nachrhPNHL8QuJJuZtQZ+SZgAOBS4yMyOiFtVvSYRaky6ncA17n4kcBIwPok/T3ffBpyW+f0+GjjN\nzE6pb9/YwQ71T3JKoluBb8cuYm/cfWPW3U7AO7Fq2Rt3f9bd92Tuzgb6xqynIe6+2N2XxK6jHiOA\nZe6+0t13AvcTJgYmirvPAN6NXce+uPvb7j4/c3sTsIgwHydx3H1L5uZ+QGtgfX37JSHY65vklChm\nNg6odPdXY9eyL2Z2g5m9CVwC3BS7nhxcBvwldhElpg+wKuu+JgEWSGbo9nGEBkfimFkrM5sPVAHT\n3b2ivv2aMiom10KeBXrV89R3CZOc/jdz/4eESU6XN3dN9dlHndcB2YuZRfuWsZc6v+PuU9z9u8B3\nM/MJfgZcWtQCM/ZVZ2af7wI73P2+ohaXJZc6E0gjHpqBmXUCHgK+nmm5J07mm+6xmfNSz5hZmbuX\n192v2YPd3T+Ry35mdgcQ7RepoTrN7CjCDNtXzAxCt8HfzWyEu68pYolA7j9P4D4itoT3VaeZfQn4\nFHBGUQpqQB4/zyRZDWSfGO/HB5fxkDyZWVvgYeAed380dj374u4bzOxJYDhQXvf52KNiemfdzZ7k\nlBjuvsDde7r7QHcfSPgFOj5GqO+LmQ3KujsOmBerlr3JjH76FjAuc0KoFCTpXNBcYJCZDTCz/Qir\nqT4euaaSZaHFdidQ4e63xa6nIWbWvaa72sz2JwzmqPd3POo4djP7A2GtmfcnObl7VbSCcmBmy4Hh\n7l7vSYuYzOwhYAiwG3gd+GpCP4CWEk7+1PwMZ7n7VRFLqpeZfRb4BdAd2ADMc/ez9/6q4shc/+A2\nwgm0O929waFvsZjZZGA00A1YA0xw90lxq/qwzMiSF4BXqe3mui5po/TMbBhwN6FB3gr4o7v/pN59\nNUFJRCRdkjAqRkRECkjBLiKSMgp2EZGUUbCLiKSMgl1EJGUU7CIiKaNgFxFJGQW7iEjK/D+UMwo2\nQcnZWQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xb0e82b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def f(x):\n",
    "    return 4*x**3 + (x-2)**2 + x**4\n",
    "\n",
    "fig, ax  = subplots()\n",
    "x = linspace(-5, 3, 100)\n",
    "ax.plot(x, f(x));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 2.804988\n",
      "         Iterations: 4\n",
      "         Function evaluations: 24\n",
      "         Gradient evaluations: 8\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ 0.46961745])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_min = optimize.fmin_bfgs(f, -0.5)\n",
    "x_min"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEACAYAAAC57G0KAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYFdW19/HvApwQBBWuGHEgiorGCQUxmNAKKJrEMTdK\nxAGjMTjGmDjmSpt7jXGKRknUOBBnVERF4wRix+FVQEVQEecRHJhEQFGG9f6xq+HQdvc5fbrq1DnV\nv8/z1MMZqk7tortX79619trm7oiISPa0SrsBIiKSDAV4EZGMUoAXEckoBXgRkYxSgBcRySgFeBGR\njMob4M1skJnNMLO3zOysRvbrZWbLzOzQph4rIiLxazTAm1lrYAQwCNgOGGxmPRrY72Lg0aYeKyIi\nycjXg+8NvO3u77v7UmAUcGA9+50CjAZmF3GsiIgkIF+A3wT4KOf5x9FrK5nZJoTAfU30Uu3U2LzH\niohIcvIF+ELqGFwJnO2h5oFFW6HHiohIQtrkeX8msGnO800JPfFcuwKjzAygE7CfmS0t8FjMTL8I\nRESK4O6Wb4cGN8IvgHeALYA1gZeBHo3sPxI4pCnHhiZk1/Dhw9NuQqIq+fref9/9gAPcN9rI/cIL\n3efN++4+ude3ZIn7jTe6b7ONe69e7s8/X7q2JqGSv3aFyPr1RbGz0Rje6BCNuy8DTgYeA6YDd7n7\n62Z2gpmdUMyxjf62ESkBd7jpJthtN+jTB957D849F9Zfv/Hj1loLjj0Wpk+H3/4WDjgAzjsPvvmm\nNO0WaapCJjp5zrYCwN2vc/frAMzsQDObamZTgB2BL3KOvQb4BlgIHBxnw0WKsWgRHHwwjBgBEybA\nOefAOus07TNatYJf/hKmToVXX4Xdd4cPPkimvSLNEUce/Hh338nddwGOAf6Z854DVe6+i7v3jq/Z\nlaOqqirtJiSqkq7vk0+gXz/o1Amefx522CH/MY1dX5cucP/9cPTR0LcvvPRSfG0thUr62hUj69dX\nCPNGFvwwsz2A4e4+KHp+NoC7/6WR/a9w9z7R8/eA3dx9biPn8MbaIBKH6dNh//3h+OPDcIw1fmuq\nycaMgRNOgJtvDucRSZqZ5b3J2uw8+OhEB5nZ68AjwKk5bzkw3sxeMLPjC2u2SLxmzIABA+BPfwpj\n5nEHd4BDDoGxY2HoUHj44fg/X6QYceTB4+73u3sP4GfArTlv9Y2GbvYDTjKzHxXXTJHivPMODBwI\nF10ERx2V7Ln22AMeeCAM2TzxRLLnEilEHHnwK7n702bWxsw2dPe57v5J9PpsM7uPUL7g6brHVVdX\nr3xcVVWlsTOJxUcfhZ77eeeFoFsKffrA6NHw85+H8fm+fUtzXsm+mpoaampqmnRMvjH4NsAbQH9g\nFjAJGJyb7mhmWwLvurubWU/gHnff0szaAq3dfaGZrQs8Dlzg7o/XOYfG4CV2CxeG4DpkCJx5ZunP\n/+ijcMwx8MwzsNVWpT+/ZF8hY/CN9uDdfZmZ1eaytwZurM2Dj96/DjgUOCqavboIODw6vAswJprh\n2ga4vW5wF0nC8uUweHDoTf/hD+m0YdAgqK6Gn/4Unnsuf469SBIa7cGXpAHqwUvMTj8dpk0Lveg1\n1lBbJJviyKLJu2hH7kQnM3vRzPYu9FiRuI0cGbJYRo8uj4B62WVhItXpp6fdEmmJ8o3BtyaMwQ8g\n3HCdzHfH4Nd198XR4x2A+9x9q0KOjY5RD15iMWUK7LMP/Oc/sN12abdmlQULQlmE6mo44oi0WyNZ\nEUcPPu+iHbXBPdIOmFPosSJxmT8/ZK6MGFFewR2gQwe4995Qv+bVV9NujbQkSU500oIfUhIrVoQc\n95/9DA47LO3W1G/HHeHyy+HQQ+HLL9NujbQU+fLgC57oBNwfTWS61cy2bUojlAcvzXHFFTBnTugl\nl7Ojjgppk8OGwW23JTOjVrIriTz4PkB1Ti2ac4AV7n5xI8e8Qxie6V7IsRqDl+aYPBl+8hOYNAm2\n2CLt1uT31VfQuzf8/vchT16kWHGMwb8AdDezLcxsTeAwYGydk2xpUbJ7NNGJqLhY3mNFmuPLL0O+\n+9//XhnBHaBtWxg1KgT4GTPSbo1kXWITnRo6NrlLkZZm2DDo3x/++7/TbknT/OAHcOGFcPjhoWzx\n2mun3SLJqmYv+EFIgVxOCOKtCAt81NKCH5KI228PaZFXXJF2S4rz61/DlluGOjkiSYkjD34PYLq7\nLzCzQYRx99x68Lu6+7xGzqExeGmS99+HXr1g3DjYeee0W1O8uXNhp53C5KyBA9NujVSaUuXBP+fu\nC6KnE4GuddvRhDaLNGrZslBA7KyzKju4A2y4IfzrX6GG/NwGl8QRKV4sefA5fgXkLnegBT8kVhdf\nHBa//t3v0m5JPAYMCGPxxx8fFgMXiVMsefAAZrYXcCyQWwG7r7t/YmadgXFmNsPdVQ9eijJ5Mlx1\nFbz4Ylj4OisuvDAs3F3bmxepT2p58Ga2IzAGGOTubzfwWcOBRe5+eZ3XNQYveS1eDD17wv/9X+Vl\nzRTi1Vdhr71CVs2WW6bdGqkEpcqD34wQ3IfkBncza2tm7aPH6wL7AK80/TJEQt54nz7ZDO4QUif/\n+Ec48shwn0EkDo0GeHdfBtTmsk8H7qrNg6/NhQfOB9YHrolKBk+KXu8CPG1mLxNuvj6kBT+kGA89\nFOqpX3112i1J1imnQLt28Oc/p90SyYo48uCfJNx8bQ18BRwX7fMucDawNrBW7bEiTfHpp+EG5C23\nwHrrpd2aZLVqFcbh//GPMFQj0lyJ5cGrHrw014oVoc7MbrvB//5v2q0pnfvuC0NSU6Zk/5eaFC/t\nPHjVg5dmGTEi1Hk///y0W1JaBx8cSjCcemr+fUUak2QevOrBS9FeeSX02m+/vTyW3iu1K64Ii3WP\nGpV2S6SSJZkHX/CxyoOXXIsXh4U7Lrus5aYMrrsu3HknDBoUcuS7dUu7RZK2ssqDb8KxGoOX1Rx/\nPCxZEm6stvRFMa68MgT6Z55pmX/JSMNSzYMv5FiRukaNCotm/+MfCu4Ap50GnTuHHHmRpoqjHnxu\nHjzAUnfvrXrw0lTvvBNywR99FNq3T7s15cEsVJvcZZcw03XQoLRbJJUkjjz4y4A3gR7Abe7eO+dY\n1YOXgnz9dViQevhw2HXXtFtTXjp3hjvuCEv8ffhh2q2RShJHHnxnYHPgIGB+bq0Z1YOXQv3qV2G9\n0jvu0NBMQy69FEaPhqeeChU1pWUrVR78bHd/AVjaUDsKbbC0TDfdFFICr79ewb0xv/89fO974V+R\nQsSdB1+X6sFLoyZPDot33HtvqMMiDTMLpQwefTRkGInkE1sefAMKqgcvLdMnn8Ahh4See48eabem\nMnToAPffH264br11qLAp0pB8AX4msGnO800JvfiCuPsn0b+zzew+wpCPFvwQliwJU/JPOAEOOijt\n1lSW7beHG28MN6UnTYJNND+8RUhiolMbwk3W/sAsYBL1FAyL9q0GFtbeZDWztkBrd18Y1YN/HLig\nbslg3WRtedzh6KNDkL/rLo27F+uii2DMGKipCTNfpWUp5CZrowE++pD9gCtZlct+UW4evJl1IWTX\nrEdIo1wIbAf8F2ECFIS/FG5394vq+XwF+Bbmj3+EceNgwgQFpuZwD6mT8+aFCpRt8v09LpkSRxYN\n5M+D70gYtlkL+D9338zdF6kevNTnuutCr/3BBxXcm8ss3L/45hs46SQt2i3flVgevOrBS11jx4Yx\n96efhq22Srs12bFwIfz4x+GG9f/8T9qtkVJJOw9e9eBlpcceg+OOC0FewT1e7dvDww+H1Mkrrki7\nNVJO8o3a1ZcHv3uBn92cYyVDJkyAIUNCel+vXmm3Jps23hieeAL69YM11wxDNiJJ5sFr3EX4z39C\nbfd77oG+ffPvL8XbbLPwy7RfP2jdGn7zm7RbJGlLMg++4GOVB59NY8eGGjN33QX6kpZGt24hyA8c\nCF98EWYJKw01G8otD76gY3WTNZtuuQXOPDNky2hYpvRmzoR994X99oNLLlGQz6JU8+DdfVF9x9bz\n+QrwGbJiBVx4IdxwQ7ixuu22abeo5Zo3D37yk9Crv+EGaNs27RZJnGIJ8ElTgM+OxYvDDNWZM8MM\ny403TrtF8vXXIXvpjTfCTe6uXdNukcQllolOZjbIzGaY2VtmdlYD+1wVvT/VzHbJef19M5tmZlPM\nbFLTL0EqxRtvhJuo7duHqfMK7uVhnXXgttvgF7+A3r3hySfTbpGUUqMBPpqsNAIYRCg/MNjMetTZ\nZ39gK3fvDvyasIpTLQeq3H2XOis9SUa4hz//99wzTGK66SYtRlFuzML9kH/9K6Srnn02fPtt2q2S\nUmj2RCfgAOBmAHefCHQ0s41y3tftnYyaOTPMnhwxIqRDDhumm3nlbJ994OWX4bXX4Ic/hKlT026R\nJC2OBT8a20cLfmTQsmVw1VWw887wgx/AxImw3XZpt0oK0blzSF8dNiwE/DPOCKUOJJvimujUUL9t\nT3eflW/BD+XBVwZ3+Pe/4bzzYMMNQ00ZZclUHrMwP+GAA8LQTY8eoYbNscfCGmuk3TppSBJ58H2A\nancfFD0/B1jh7hfn7HMtUOPuo6LnM4B+7v5Znc8aDizKXZQ7el1ZNGXOPZT3veACWLAgpEEecICG\nY7Ji8mQ491x47z04/3w4/PBQ7kDKWxxZNC8A3c1sCzNbEzgMGFtnn7HAUdEJ+wBfuPtnZtbWzNpH\nr68L7AO8UsR1SEq++iqsHLTDDuFP+WHDwrjtgQcquGdJr17hF/g//wk33xzy5v/8Z5gzJ+2WSXM1\nOkTj7svM7GTgMVZNVno9d6KTuz9sZvub2dvAYmBodHgXYIyFSFC74Mfj3z2LlJPly8PQy623hlz2\nvn3hyiuhf38F9azbe++wTZ0avuZbbRVKTBx1FOy/P6y9dtotlKaKY8EPotc9Z3+04EfQ1DGzNMyZ\nA6NHh9WBunSB3/423DSdPh0eeggGDGg4uFfC9TVHlq+voWvbaScYORI+/DD8tXb11bDRRmHt3Btv\nhPffL2kzi5blr12hEsuDL+TYlqDcvsmWLIEXXwwrAf3mN2EB5y23DPnru+0GL7wQUunOOKOwyUrl\ndn1xy/L15bu29daDoUPD5Kh334Wf/zwM5ey+O2y+ORx5JPztb+Evvi+/LE2bmyLLX7tC5cuiWZkH\nD2BmtXnwuQXDVsuDN7OOUX2abgUcKzH7+mv49NOwzZoFH3wQelzvvgszZsDHH0P37rDrrtCzZ5jG\nvvPOWs9TGrfhhmGS1JAh4ab7m2+GwP7SS3DnnfDKK7D++iGraqutwi+AzTcPpRG6dAmdhXbtNMxX\nanEs+NFQHvz3Cji2pGbNClPqi1Ffok/ua7WP3VdtAG+/DY88suq9FStWbbXPly8PueXLl696vHRp\n2L79NmxLloTt669DzZfFi2HRotBzWrAA5s+HuXPD8V26rPqh2mKLcNOsf/+QDtetm1LhpHnMYJtt\nwlZrxQr46KPw8/XWW6FjMWVK+Jmr7XAsXQobbBC2Dh3CXwjt24e1edddNxRDW3vtMBN6rbVCJs8a\na4StTZuwtW69amvVavXNbPXt7bdDwbva57Vtz/0l09AvnGJ+EW28cRmmDbt7gxtwKHB9zvMhwNV1\n9nkQ6JvzfDywayHHRq+7Nm3atGlr+tZY/Hb3WBb8qLtP12ifNQo4Nm8ep4iIFCexPPgCjxURkYQk\nlgff0LFJXoyIiKyS+oIfIiKSjEImOomISAVSgBcRySgFeBGRjFKAFxHJqKIDvJndZGafmVmDJYAb\nWoxbRESS15we/EhCIbF65VmMW0REElZ0gI+W3pvfyC75FuMWEZEEJTkGX18Rsq4Jnk9ERHIkXSS2\nbp2Z78yqMjPNtBIRKUK+Wl5JBvj6ipDNrG/HLM+mra6uprq6Ou1mJEbXVznmzw9rArz5ZijpO3Zs\nNeutV81HH4Vyvu3ahZWbOneGTp1CWd/114eOHVeV9m3XLmxt28I666xe3re2xG+bNqtK/LZunV4N\n+Cx97epjBfzHJhngxwInA6PqFCETkYR9/jlMnBi2KVNg2rSwbsC228LWW4dFX7bcMizPuNlmoZa5\n1lzNnqIDvJndCfQDOpnZR8BwQongfItxi0jM5s4Ny+k9+WTYZs+G3r3D8nrHHx/WWd1887AwRq3q\naujXL7UmSwkUHeDdfXAB+5xc7OdnRVVVVdpNSJSuLz3vvgv33BMWRp86FaqqYO+94cQTYYcdVg/m\n9Snna4tD1q+vEKlXkzQzT7sNIpVi7ly47Ta44w547z049FA48MAQ3DXE0rKYWd6brArwImXOHZ55\nBq69Fv79b/jpT+HII8M6u1osveVSgBepYMuWwb33wuWXwxdfwEknhcC+wQZpt0zKQSEBXr//RcrM\nsmVhCOZPfwrZLX/8Y+i15xtTF6mrWd8yZjbIzGZEBcXOquf9Tmb2qJm9bGavmtkxzTmfSJa5w+jR\nsP32cMMNYXv6aTjgAAV3KU7RQzRm1hp4AxhAmMA0GRicu+6qmVUDa7n7OWbWKdp/I3dflrOPhmik\nxZs0CX73O1i0CC69FAYMSG+CkFSGQoZomtMv6A287e7vu/tSYBRwYJ19PgHWix6vB8zNDe4iLd28\neSFP/aCD4Nhj4cUXYeBABXeJR3MCfH3FxDaps8/1wPZmNguYCpzWjPOJZIZ7SHfcbruQ3vj66yHA\nt26ddsskS5pzk7WQcZVzgZfdvcrMtgTGmdlO7r6wGecVqWiffAK//jV8+CGMHRtmnIokoTkBvm4x\nsU0JvfhcPwQuBHD3d8zsPWAb4IXcnXILAlVVVWkGmmTW3XfDKaeEAH/vvaE4l0ghampqqKmpadIx\nzbnJ2oZw07Q/MAuYxHdvsv4VWODuF0SLfbwI7Oju83L20U1WybzFi+HUU8OEpdtug1690m6RVLpE\nb7JGN0tPBh4DpgN3ufvrZnaCmZ0Q7fZnYDczmwqMB87MDe4iLcG0abDbbrB8ebiJquAupaKZrCIJ\nuuUWOOMMuOIKGDIk7dZIlmgmq0hKvvkm5LWPGwc1NWHykkipKcCLxOzzz+Hgg8PKSJMnQ4cOabdI\nWipNgBaJ0dSpIe1xwAAYM0bBXdKlHrxITMaOheOOgxEj4Be/SLs1Is3owecrNBbtU2VmU6JCYzVF\nt1KkzF11FQwbFuq1K7hLuSgqi6bAQmMdgWeBfd39YzPr5O5z6vksZdFIxVq+PNxMHT8+BPcttki7\nRdJSJJlFs7LQWHSi2kJjr+fs80vgXnf/GKC+4C5SyZYsgSOOgPnz4dlnoWPHtFsksrpih2gKKTTW\nHdjAzJ40sxfM7MgizyVSdubPh333hTXWgEceUXCX8lRsgC9kTGUNoCewP7Av8D9m1r3I84mUjZkz\n4cc/hp49w8pLa62VdotE6lfsEE0hhcY+Aua4+9fA12b2FLAT8FbdD1OxMakUb70F++wDJ5wAZ52l\nuu1SOiUrNlZgobFtgRGE3vtawETgMHefXuezdJNVKsKUKfCTn4S1Uo87Lu3WSEuX2E1Wd19mZrWF\nxloDN9YWGovev87dZ5jZo8A0YAVwfd3gLlIpnnoKfv5zuPZaOOSQtFsjUhgVGxPJ4+GH4Zhjwnj7\ngAFpt0YkSHpNVpHMGzUKhg4Ns1QV3KXSKMCLNOC660Kp3/HjoU+ftFsj0nSqRSNSj0sugWuugf/8\nB7baKu3WiBRHAV4khzucdx7cdx88/TR07Zp2i0SKl2ixsWi/Xma2zMyUeyBlbcUKOPFEePxxBXfJ\nhqJ68FGxsRHkFBszs7G5efA5+10MPApoSoiUraVL4eijYdYsmDAB1lsv7RaJNF+xPfiVxcbcfSlQ\nW2ysrlOA0cDsIs8jkrivvoIDD4SFC0NdGQV3yYrEio2Z2SaEoH9N9JKS3aXszJ8PAweG5fXGjIF1\n1km7RSLxSbLY2JXA2dEsJkNDNFJmaouG9ekDI0eGypAiWZJksbFdgVEWqjF1AvYzs6XuPrbuh6nY\nmJTa66/DfvuFm6p/+IOKhkn5K6tiY3X2Hwk86O5j6nlPpQqkpJ59NtSTuewyOFKrFEiFSrXYWDGf\nK5K0e+6Bk06CW28NC3aIZJmKjUmL4A6XXgpXXw0PPQQ77ZR2i0SaJ8k1WUUqxrffhl77pEnw3HOa\nwCQthwK8ZNrs2aGOe4cO8Mwz0L592i0SKR1Vk5TMmjYNdt8d+vaF++9XcJeWRz14yaTbboPTT4e/\n/Q1++cu0WyOSjmb14PMVHDOzI8xsqplNM7NnzWzH5pxPJJ9vv4VTToELLoAnnlBwl5at6B58gQXH\n3gV+7O4LzGwQ8E9ASydIIt5+GwYPhk02gcmToWPHtFskkq7m9ODzFhxz9+fcfUH0dCKg/AVJxO23\nwx57hIqQ992n4C4CzRuDr6/g2O6N7P8r4OFmnE/kO+bMCSmQ06bBuHGw885pt0ikfDSnB1/w7CQz\n2ws4FmhwYRCRpho7FnbcMeS1v/SSgrtIXc3pwRdScIzoxur1wCB3n1/fB6nYmDTFxx/DaaeFXvtd\nd8GPfpR2i0SSV7JiY1BYwTEz2wyYAAxx9+cb+ByVKpCCfPst/P3vcOGFYVjmnHNg7bXTbpVIOhIt\nVVBgwbHzgfWBa6KywUvdvXex55SWyR0eeCCU9e3ePVSD3GabtFslUv5UbEzK2pNPwvnnh5WXLr9c\nFSBFaqnYmFQkdxg/PgzFzJwJw4eH/PbWrdNumUhlUYCXsrFkSbhp+te/wooV8PvfwxFHQBt9l4oU\nRT86krrXXoMbbgj1Y3r2hEsugX320TJ6Is2lAC+pePdduPtuGDUqlPQdOjTUa+/WLe2WiWSHbrJK\nSXzzDTz/PDzySFhR6fPPw7qogwfDnntqfF2kqQq5yZpoNclon6ui96ea2S7NOV8laurEhErT0PXN\nnRuCeXU19O8PG24YxtRbtw7DMZ9+CtdeC/36lXdwz/LXL8vXBtm/vkIUHeBzqkkOArYDBptZjzr7\n7A9s5e7dgV8D1zSjrRUpy99k7vDIIzVMnBgWsT7vPPjZz2DzzcNQy6WXhp77b38bZp9OnhwyY/r0\ngVYVstRMlr9+Wb42yP71FaI5Y/Arq0kCmFltNcnccsEHADcDuPtEM+toZhu5+2fNOK8kaNkyWLgQ\nvvgibPPmhTHyOXPgs89g1iz45BP48EP44IMwu3TCBNh667ANHRrqw3z/+5UTxEWyKulqkvXt0xVI\nJcBPmRKGDZKWe0vhqadCr7Wh9xt7XPu89nHttmLFdx8vXx4eL18etmXLVt+WLg3BeOnS0KtesiRs\nX38NX30FixfDokVhn3btYP31w9axI3TuHLb/+q+wBN7GG8Nmm4We+pVXhmEYESk/zalFcyihgNjx\n0fMhwO7ufkrOPg8Cf3H3Z6Pn44Ez3f2lnH10h1VEpAhJzmQtpJpk3X26Rq8V3EARESlOc0ZJXwC6\nm9kWZrYmcBgwts4+Y4GjAMysD/CFxt9FREoj0WqS7v6wme1vZm8Di4GhsbRaRETySn2ik4iIJEOJ\nbCIiGaUALyKSUQrwIiIZpQAvIpJRsQR4M9vGzKbkbAvM7FQz28DMxpnZm2b2uJl1jON8IiKSX+xZ\nNGbWijCZqTdwCjDH3S+Jqk2u7+5nx3pCERGpVxJDNAMIRcg+IqfYWPTvQQmcT0RE6pFEgD8cuDN6\nnFs58jNgowTOJyIi9Yh1iCYqWTAT2M7dZ5vZfHdfP+f9ee6+QZ1jNNNKRKQIia7oVI/9gBfdfXb0\n/DMz6wJgZhsDn9d3kLtndhs+fHjqbdD16fpa2rW1hOsrRNwBfjCrhmcgFBs7Onp8NHB/zOcTEZEG\nxBbgzWxdwg3WMTkv/wUYaGZvAntHz0VSY2bf2USyqjn14Ffj7ouBTnVem0cI+i1WVVVV2k1IVGVe\nX+6ft40H+Mq8vsJk+dog+9dXiNhuskaTmG4Atif8BA0F3gLuAjYH3gd+4e5f1DnO42qDSD6hx756\ngM/9/quvR6/vTylHZoaX8Cbr34CH3b0HsCMwAzgbGOfuWwNPRM9FypznbCKVK5YevJl1AKa4+/fr\nvD4D6Ofutdk0Ne6+bZ191IOXxNQ/xp6vB9/w+yLlopQ9+G7AbDMbaWYvmdn10U1XTXSSMqAeubRM\ncQX4NkBP4B/u3pOwPN9qwzFRN10/YSIiJRJXFs3HwMfuPjl6Pho4B/jUzLq4+6eNTXSqrq5e+biq\nqkp3v6Vs6SaspKWmpoaampomHRNnFs1TwHHu/qaZVQNto7fmuvvFZnY20NHrVJPUGLwkqb4x9eaM\nwWuMXspFIWPwcQb4nQhpkmsC7xDSJFsDdwOboTRJSYECvGRVSQN8sRTgJUkK8JJVpc6DFxGRMhJb\nqQIzex/4ElgOLHX33ma2AXlmsoqISDLi7ME7UOXuu7h77+g1zWQVEUlJ3EM0dceDtGSfiEhK4u7B\njzezF8zs+Og1zWSVklIpYJFVYhuDB/q6+ydm1hkYF9WhWcndvaHl+TTRSYpVSK2ZZM4hUlqpTnRa\n7UPNhgOLgOMJ4/K1M1mfVLExiVMhaZDNTZNsyvH1t291+n6XOJQsTdLM2ppZ++jxusA+wCtoyT4R\nVOxM0hLXEM1GwH1Rb6UNcLu7P25mLwB3m9mviNIkYzqfSGw0BCNZpZmsUtHiGKKJc4insPatTt//\nUoySz2Q1s9ZmNsXMHoyeb2Bm48zsTTN7PFrWT6Ro2ciS0ZCNlEbcefCnAdNZ9Z2riU6SAAVIkULE\nFuDNrCuwP6GiZG3XShOdRERSEmcP/grgD8CKnNc00UlEJCWxZNGY2U+Bz919iplV1bePJjqJiBQv\ntYlOZvZn4EhgGbA2sB4wBuiFJjpJjMohayb/RKm6is/CEWlIybJo3P1cd9/U3bsBhwMT3P1INNFJ\nWiTdBJbykNSCH7Xf2X8BBprZm8De0XMRESkBTXSSilIZQzTxTZQSaUgpa9GsbWYTzexlM5tuZhdF\nr2uik0ge2Zi8JeUorjH4JcBe7r4zsCOwl5ntiSY6SQvQ/OCsMXtJRmxj8O7+VfRwTaA1MB9NdJIW\nId4ArR7+AnP0AAAGsElEQVS9xCXOmaytzOxlwoSmJ939NTTRSaQI6tFLPGJb0cndVwA7m1kH4DEz\n26vO+w1OdBJpiHqwIsWLc8k+ANx9gZn9G9gV+MzMuuRMdPq8vmM0k1UaF+8SfCKVKM2ZrJ2AZe7+\nhZmtAzwGXADsC8x194vN7Gygo7ufXedYpUlKg8oxLbKc6s1Ly1VImmRcPfiNgZvNrBVhXP9Wd3/C\nzKagFZ2kCTQkIxIfTXSSsqIeu3rwUpiSr+gkIiLlI66ZrJua2ZNm9pqZvWpmp0avayarSIzq5shr\nSEsaE1cPfilwurtvD/QBTjKzHmgmq0gClCcvhYmrVMGn7v5y9HgR8DqwCZrJKiKSmtjH4M1sC2AX\nYCKaySoikppYJzqZWTvgXuA0d1+YOz6oJftapvrGiJUVUjr6/8+O1CY6AZjZGsBDwCPufmX02gy0\nZF+LVl/aY2Nfb6VJNr3efFP/P/Xzlg2lrAdvwI3A9NrgHtGSfdIoZYSIJCeuUgV7Ak8B01jVXTgH\nmATcDWxGNJPV3b+oc6x68BmWrwepHrt68FKcQnrwmskqiVKAV4CXZJRyiOYmM/vMzF7JeU2TnERi\nEPcQlobFWo640iRHAoPqvKZJTiKxaHxiU9MDtiZKtRRxTXR6mrBEXy5NchIpCQVsqV/sC37k0CQn\nkRRo2EVqJRngV8q3XJ8mOrUsCkBJq3vTVrIg7YlOWwAPuvsO0fO8k5yi/ZRFkyH1B+/yykrRc2XV\nZEHa9eA1yanF0phwJVFWTXbFNdHpTqAf0Ikw3n4+8AB5JjlFx6oHnyHKa8/C89Xp57M8aaKTJE5D\nMtl/rp/P8pT2EE1tIwaZ2Qwze8vMzkr6fJIGDcmIlKNEA7yZtQZGECZBbQcMjlZ6ajGaetc7bvmW\neCvmfY3XZkVN2g1IVNo/e+Ug6R58b+Btd3/f3ZcCo4ADEz5nWSmPb7J8PezV3/9uAK97vHrs2VBT\n0F7N6SCkqTx+9tKVdB78JsBHOc8/BnZP+Jxla+TIm3nggfErn7dqBSeeeDQDBgwoaTvy/xAqj1py\n5ft+0PdLuUo6wKuLl+P//b/JPPDAbau9NmDAHrEG+MJrkaw8IrZzi0D+78HvVsds+H1pnkSzaMys\nD1Dt7oOi5+cAK9z94px99NUUESlCqmmSZtYGeAPoD8wiLAAy2N1fT+ykIiICJDxE4+7LzOxk4DGg\nNXCjgruISGmkPtFJRESSkfhEp0KZ2Slm9rqZvWpmF+c/ovKY2RlmtsLMNki7LXExs0ujr9tUMxtj\nZh3SblMcsjxBz8w2NbMnzey16Oft1LTblAQza21mU8zswbTbEicz62hmo6Ofu+nRvc56lUWAN7O9\nCAuE7OjuPwAuS7lJsTOzTYGBwAdptyVmjwPbu/tOwJuExdYrWguYoLcUON3dtwf6ACdl7PpqnQZM\nJ3vZfH8DHnb3HsCOQIPD3mUR4IFhwEXRZCjcfXbK7UnCX4Ez025E3Nx9nLuviJ5OBLqm2Z6YZHqC\nnrt/6u4vR48XEQLE99JtVbzMrCuwP3ADGcoFjv5C/pG73wThPqe7L2ho/3IJ8N2BH5vZ82ZWY2a7\npd2gOJnZgcDH7j4t7bYk7Fjg4bQbEYP6JuhtklJbEhWt47AL4ZdzllwB/AFYkW/HCtMNmG1mI83s\nJTO73szaNrRzSVZ0AjCzcUCXet46L2rH+u7ex8x6EcoMf79UbYtDnus7B9gnd/eSNComjVzbue7+\nYLTPecC37n5HSRuXjKz9SV8vM2sHjAZOi3rymWBmPwU+d/cpZlaVdnti1gboCZzs7pPN7ErgbEKJ\n9np3Lgl3H9jQe2Y2DBgT7Tc5uhG5obvPLVX7mquh6zOzHxB+606NZu11BV40s97u/nkJm1i0xr52\nAGZ2DOHP4f4laVDyZgKb5jzflNCLzwwzWwO4F7jN3bO2GM8PgQPMbH9gbWA9M7vF3Y9KuV1x+Jgw\nGjA5ej6aEODrVS5DNPcDewOY2dbAmpUU3Bvj7q+6+0bu3s3duxG+QD0rJbjnY2aDCH8KH+juS9Ju\nT0xeALqb2RZmtiZwGGGFskyw0NO4EZju7lem3Z64ufu57r5p9PN2ODAhI8Edd/8U+CiKkwADgNca\n2r9kPfg8bgJuMrNXgG+BTHwxGpC1P/+vBtYExkV/oTzn7iem26TmaQET9PoCQ4BpZjYleu0cd380\nxTYlKWs/c6cAt0edj3eAoQ3tqIlOIiIZVS5DNCIiEjMFeBGRjFKAFxHJKAV4EZGMUoAXEckoBXgR\nkYxSgBcRySgFeBGRjPr/SYBJHv8M/N8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x15241390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# create a (continous) random variable with normal distribution\n",
    "Y = stats.norm()\n",
    "\n",
    "x = linspace(-5,5,100)\n",
    "\n",
    "fig, axes = subplots(3,1, sharex=True)\n",
    "\n",
    "# plot the probability distribution function (PDF)\n",
    "axes[0].plot(x, Y.pdf(x))\n",
    "\n",
    "# plot the commulative distributin function (CDF)\n",
    "axes[1].plot(x, Y.cdf(x));\n",
    "\n",
    "# plot histogram of 1000 random realizations of the stochastic variable Y\n",
    "axes[2].hist(Y.rvs(size=1000), bins=50);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 1.0, 1.0)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y.mean(), Y.std(), Y.var()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1.193733722708372, 0.23272385793970249)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# t-test example\n",
    "t_statistic, p_value = stats.ttest_ind(Y.rvs(size=1000), Y.rvs(size=1000))\n",
    "t_statistic, p_value"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
